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AN EFFICIENT CODE FOR THE SIMULATION OF 


NONHYDROSTATIC STRATIFIED FLOW OVER OBSTACLES 

Gregory G. Pihos and Morton G, Wurtele 
Department of Atmospheric Sciences 
University of California, Los Angeles 

INTRODUCTION 

Background 

The gravity wave is the subject of a voluminous literature containing 
theoretical and/or observational studies. However, the numerical simula- 
tion of this phenomenon has received much less attention than has climate 
modeling or numerical weather prediction. A two-dimensional , nonlinear, 
nonhydrostatic model (Foldvik and Wurtele, 1967) produced realistic 
results, but was too expensive for operational use. Various linear models 
(Danielsen and Bleck, 1 970; Vergeiner, 1971) avoided this difficulty, but 
could be considered reliable only when reproducing wave-like features and 
not when simulating turbulence-generating, wave-breaking patterns. Some 
highly successful computations are those of Klemp and Lilly (1978), which 
are applicable when the disturbance generated satisfies the quasi-hydro- 
static approximation. 

When the prediction of areas of clear-air turbulence (CAT) is the 
chief emphasis of a study, it is essential to retain both nonlinear and 
nonhydrostatic effects in any numerical model. Since 1966, computers and 



computational techniques have been developed to an extent permitting the 
formulation of a model like that of Foldvik and Wurtele, but efficient 
enough to run at low cost. As a consequence, such a model can be used 
(1) to study sensitivity of results to input data; (2) to test the im- 
plications of a great variety of idealized initial and boundary conditions; 
and (3) to simulate easily and cheaply in real time the gravity-wave and 
CAT patterns associated with any operationally analyzed or predicted syn- 
optic situation. The model and the code of the present study have been 
developed with all three of these purposes in mind. 

Gravity Waves and CAT 

Before describing the model in detail, it may be advisable to clarify 
some ideas and concepts involved in the relation between gravity waves and 
clear air turbulence. Although a number of subtle dynamic considerations 
are involved in the stability of stratified shear flow, for the purposes 
of this study we shall proceed from the assumption that a Richardson number 
of 0.25 is the marginally critical value for the stability of an incom- 
pressible Boussinesq flow to small perturbations. Normally the initial 
conditions assumed for the model will be characterized by Richardson num- 
bers many times larger than the critical. By one means or another-- 
represented in the model by flow over an obstacl e--thi s flow is disturbed, 
and disturbed flow will contain areas in which the Richardson number is 
reduced from its initial value and areas in which it is increased. The 
dynamic/kinematic mechanism of this Richardson number modification-by- 
deformation is subject to various semi -quantitative explanations. Some 
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interpretations are reviewed by Pao and Goldburg (1969). The most widely 
accepted explanation of CAT generation is that large amplitude gravity 
waves resulting from flow over mountain ranges can and do generate local 
regions in which the Richardson number falls below the critical value, and 
moderate to severe CAT results. 

Steady state linear or nonlinear solutions have been constructed for 
a number of highly idealized conditions. However, mathematical analysis 
cannot predict with any precision whether and where areas of subcritical 
Richardson number will occur from an arbitrary disturbance in an arbitrary 
flow field. Thus, for the purposes of the present work, we rely upon the 
numerical model exclusively to make these predictions. No attempt is made 
here to verify any particular theoretical interpretations of CAT formation 
this must be the goal of further study. It should also be emphasized that 
qualitative features, such as wavelength, rotor formation, trapping, and 
upward propagation, are v/ell simulated by various gravity wave models in 
the literature, but that there exist few quantitative comparisons with 
observational data. Further comparison to actual measurements will be 
required to measure the reliability of this model, or any other, to 
simulate nature. To this end, use of this model is welcome, and we have 
endeavored to make the code as understandable and as versatile as possible 
In addition to the description of the model in the following section, a 
documentation of the code is given in the appendix. 


THE TWO-DIMENSIONAL, BOUSSINESQ MODEL 

The model must include buoyancy as the primary restoring force for 
any disturbance of the free stream. Hov/ever, dynamic compressibility — 
the effect resulting in acoustic waves — is not significant in the study 
of CAT. Thus, incompressibility is assumed, but in a manner that retains 
the static effect of compressibil ity. Thus temperature, potential 
temperature, density, and pressure in the undisturbed atmosphere may be 
realistically represented in the initial conditions of the model. Further, 
the Boussinesq assumption is made, neglecting the kinematic effect of den- 
sity variation (that is, where density multiplies velocity) while retaining 
the full dynamic effect of buoyancy (where density multiplies gravity). 

The two-dimensional, Boussinesq model greatly facilitates computational 
ease and speed. Only tv/o variables, vorticity and density, are directly 
advanced in time. The third variable, the streamfunction , is obtained 
at each time step by solving a Poisson equation (eq. (7c) below). The 
method of solution consists of applying a fast Fourier transform in the 
horizontal, then utilizing a one-dimensional marching solution in the 
vertical. This noniterative procedure is at least an order of magnitude 
faster than iterative relaxation techniques. Another advantage of this 
model is the existence of an energy integral for arbitrary mean density 
profiles, such as upstream inversions, (in contrast, non-Boussinesq models 
present computational difficulties when the stability profile is varying 
rapidly.) These factors will permit the program to be used frequently with 
sounding data, or in theoretical profiles. A description of the input/output 
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options is given in the appendix. 


Equations Solved by the Model 


Vie begin with the following set of equations: 
The equation of motion: 


dv -1 

§+fkxv= ivp + g 


(U) 


The equation of state: 

p=pRT (lb) 

The equation of continuity: 

^+pV-v=0 (Ic) 

and the thermodynamic equation: 

f = -R|f(1np) (Id) 

Molecular viscosity and thermal conductivity will be considered to be 
unimportant. 

This set of four equations in four unknowns may be simplified by 
making several assumptions. First, we will concentrate on mountain- 
induced gravity waves, and will limit our consideration to length scales 
of motion which are small compared to cyclone scale motions, so the 
effect of rotation will be ignored. Second, the time scale of the motion 
is small compared to the time scale of radiative heating, so the adiabatic 
assumption, = 0, may be made. The system of equations then becomes: 


-iv 

dt p '' 


P + 


->• 

g 


(2a) 


5 


p=pRT 


(2b) 


^ + pV*v=0 (2c) 

% W '1" T) - (in P)=0 (2d) 

A further simplification results from assuming the fluid density to be 
incompressible but not homogeneous. This permits us to allow density grad- 
ients in the vertical and buoyant restoring force without the unnecessary 
computation of acoustic motions. Both terms in equation (2c) then become 
equal to zero. This permits us to express the system of equations (2) as a 
system of three equations in three unknowns: 


dt - p ^ p + g 

(3a) 

dt ^ 

(3b) 

V*v=0 

(3c) 


The Boussinesq approximation states that the kinematic effects of density 
gradients are negligible compared to the buoyancy effects of density 
gradients in the equation of motion. The nonlinear equation (3a) then 
becomes : 


dv 


Po S 


pg 


(4) 
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v/here is an undisturbed reference density. 

Although significant airflow is deflected around long mountain ranges 
such as the Sierra Nevada, even when the wind is normal to the ridgeline 
(Holmboe and Klieforth, 1957), we will concentrate on the flow passing 
directly over the mountain. For this purpose, the cross section of the 
range is sufficiently uniform such that the flow may be assumed to be 
two dimensional. Taking the curl of equation (4) with the operator (j^'Vx) 
in a left-handed x-z coordinate system yields the vorticity equation; 


||=-V.(,v) 


1 ^ 

Po 


Similarly, = -V«(pv) 


(5a) 

(5b) 


v/here i, = l*Vxv = 


3u 

3z’ 


and u and w are the horizontal and vertical 


components of v. 

The streamfunction Tp for two-dimensional incompressible flow may be 

2 2 

defined by |^ = w and = u. Then, C = where 
ax az 

and the system of equations (5) becomes: 


_ 8^ 84) 8^ dtp ^ 8p 

8t ~ 8x 8z ” 8z 8x " p 8x 


(6a) 


8p ^ 8p 8 t|/ _ 8p dtp 

8t “ 3x 8z “ 8z 8x 


(6b) 



(6c) 
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or introducing the Jacobian operator, J(a,b) = 


3a 3 b 
3z 3x 



g 8p 
■ P„ 9x 

(7a) 

It ^ 


(7b) 



(7c) 


The lower boundary condition is that the surface is a streamline: 

i|;(x,h(x)) = constant (7d) 

where h(x) is the height of the barrier. The upper boundary condition is 
that the kinetic energy becomes vanishingly small with elevation: 

^ p{u^ + w^) = Y 0 as 2 “ . (7e) 

The system of equations (7) constitutes the model on which the program 
is based. The scheme for their numerical solution and the associated code 
are fully described in the appendix. The assumptions made are justified 
and discussed below. In subsequent sections, numerical simulations with 
the model are compared with selected special cases for which analytical 
solutions can be obtained. 


A Justification of the Incompressible and Boussinesq Approximations 


The effect of these two assumptions may be seen by examining the ver- 
tical velocity, profil e of the steady state, linearized perturbation equa- 
tions. Defining 0 = T(— and IT e c {-^ where k = R/c , the compres- 

a Vp / p> 

sible system of equations (2) may be reexpressed as: 
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(8a) 



-6Vn + g 




(8b) 


de 

dt 


= 0 


(8c) 


2 

v/here = CpK0n/c^ = CpRT/c^. In a two-dimensional steady-state system 
(where fr-= 0), these become: 

o t 


U 


8u 

8x 


” Sz ® to 


(9a) 




9w 

9z 


= -0 


M. 

92 


- g 


. (9b) 


0(u 


9x 


+ w 


- r 

9z^ " "^s ^9x 9z^ 


(9c) 


u 


M. 

9x 


+ V/ 


9z 


= 0 


(9d) 


Now, assume that each variable q(x,z) in the system may be expressed as 
the sum of a perturbed part q'(x,z), and an unperturbed part q(z); that 
is, q(x,z) = q'(x,z) + q(x,z). Further assume that w(z) = 0, and that 
the unperturbed state is in hydrostatic balance, that is, g = -0 • 

The linearized perturbation equations for the compressible system are then 




M - A M' 

9z “ 9x 


(10a) 


9 . 


(10b) 


- 3w'_ S II', a . M 

“ 3x ® az ® az 


- 30' , 36 a 

“ 3F' ° 


(10c) 


u 9iix _ - 2, 


(lOd) 


This system is then solved for w' : 


Aia 4-'- ly in(-) If'-^ 'M? t> '"S -Iy I"'-' i - “ 

3Z^ p u p u 9z^ 


where m = 1 - . Usually, u « c , which means m-1 , so the equation 

S 

for w' in the compressible system becomes: 


a^wV 8^w' 3w' . (9S s 9u 1 9^U\,,> _ a 

— o- + — o -Sa7 + (-7 + - - T — ~ ^ 

9x^ 9z^ u a a 9z^ 


1 9p c- - 1 99 j e 9 
where s = — ^ = T ^ s = S + — ^ . 


To analyze the effect of the incompressible assumption, v;e will derive 
a similar equation for w' from the incompressible system of equations (3). 
In a two-dimensional steady-state system, these equations have the form; 


9u , 9u -1 9p 

' Fa? 


(13a) 


9w . . 9w 

" 97 ^ 


= -li£ - 


(13b) 


10 


(13c) 


u 


9^ 

3x 


+ w 



+ 9w 
9x 8z 


= 0 


(13d) 


As before, assume that each variable q(x,z) may be expressed as 

q(x,z) = q'(x,z) + q(x,z), and assume that gp = ~|^ . Then the linearized 

o Z 

perturbation equations for the incompressible system are: 



3u _ -1 8 p ' 
3Z ~ - 3x 

(14a) 


P 



. ^ 
3z 

(14b) 

p 

P 


= 

ll= 0 

(14c) 

1^-'+ ^' = 
3x 3z 

0 

(14d) 


The equation for w' in the incompressible system is then: 


1 8^u^ 


&^w' ^ 3w‘ , , S 3u I o u\..i ^ 

— ^ 7 - S — + ^ ?)w = 0 

3x 3z 3z u^ u u 3z^ 


(15) 


The only difference between the steady-state compressible and incom- 

9 3w Q s w * 

pressible equations appears in the terms and . So atmospheric 
motions may still be modeled by the heterogeneous incompressible model 

by replacing the frequency of oscillation of the incompressible fluid 

1 /2 1 /2 
(gs) ' by the Brunt-Vai sala frequency N = (gS) ' of atmospheric 
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buoyancy oscillations. 

Now, to discuss the effect of the Boussinesq approximation, we will 
derive an equation for w' from the system of equations (4), (3b), and (3c). 
The two-dimensional steady-state equations are: 


u an + „ = -1 an 

8X " 8z Pq 8x 

(16a) 

^ 8x + 8z 8z 

(16b) 

The complete set includes equations (13c) and (13d). 
turbation equations are: 

The linearized per 

0 anV al = -I in' 

8x 8z Pq 8x 

(17a) 

- 8w' _ -1 p'g 

3X 3X 

(17b) 


The complete set of linearized equations includes (14c) and (14d), These 
result in the following equation for w' in the Boussinesq system: 


7 2 

8 w'. 8 w’ 


8X 


8Z 


+ (^- = 0 
u u 8Z 


(18) 


where . 

0 Pg 3z 

The incompressible Boussinesq and non-Boussinesq w' equations may be 
made similar to each other by transforming the non-Boussinesq equation into 
the form: 
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(19) 


2 . 2 . 

8 CO , 3 CO , 

2 2 
3X 3z 


?• 


_s 

u 


3z 


i ^ 

u 8z^ 


- Is2 
4 ^ 


3Sn , 
3F)“ = 


where w' w' . In most atmospheric profiles, the ^ and -7^—^ 

^0 g I" 12 1 3s ^ ^ 

terms dominate over the ~ 7 ^, -jr s , and ^ terms, so the dynamics of w* 

“ d Z ^ d oZ 

in the Boussinesq system are similar to the dynamics of w* in the non- 

Boussinesq system. Since w' = co'j the dynamics of the two systems 

P 

are similar to the degree that p(z) remains constant. This is true essen- 
tially for shallow atmospheric systems (Ogura and Phillips, 1962). 

To sum up, the reason for using the incompressible assumption and 
Boussinesq approximation is to simplify the system of equations to be 
solved. In order to retain the dynamics of the compressible atmosphere, 
the stability S will be used, and vertical density gradients will be re- 
tained everywhere. The total percentage variation of density in the fluid 
will be kept small, since p vn'll then have the same scale height as the 
potential temperature. Then the compressible system w' equation (12) 

applies, and may be similarly transformed using w' = (— ) ^ w' to yield: 

P 

^'+ ^'+ F(z)co' = 0 (20) 

8x Sz 


where F(z) = ^ 


I3u 1 2j^l3s 
G ^ 2 ^ ' 


Techniques for Analytical Solution of the Linear Problem 

The mountain wave problem consists of solving equation (20) with the 
appropriate boundary conditions in the upper half plane z ^ 0. The lower 
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boundary condition consists of tangential flow to the barrier. 


(x,z) = [u{z) + u’(x,z)] at z = h(x) 


where h(x) is the height of the barrier. Assuming that h is small, the 
linearized version of equation (21) is written as: 


w'(x,0) = 0(0) 


The upper boundary condition is that the kinetic energy, pw'^/2, vanish at 
2 2 1/2 

r = where r = (x + z ) ' , In terms of w'(x,z), these conditions be- 
come : 


i' (x,0) = u(0) 


dh(x) 


, and 1 im co' = 0 . 


Now, assume that o)'(x,z) may be expressed as a sum of individual wave 

oo ^ 1 k X 

components, m'(x,z) = / w(k,z)e dk, and express the ground terrain as a 

^oo 

sum of Fourier components, h(x) = /°°h(.k)e^ '^^dk. Since the system is 

— oo 

linear, the behavior of a single wave component may now be examined. The 
CO system becomes: 


d m(k,z) 


+ [F(z) - k'^]co(k,z) = 0 


(24a) 


with boundary conditions: 


)(k,0) = iku(0)h(k), and limco(k,z) = 0. 
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The general solution is oiCk.z) = c^(k)w-|(k,z) + C 2 (k)w 2 (k,z) , where c-|(k) 

/N 

and C 2 (k) are arbitrary constants to be determined by lim [ai(k,z)l = 0. 

z-x» 

Then w(k,z) = C2(k)o32(k,z) , where w^CkjZ) is a linear combination of w^(k,z) 
and co 2 (k,z), and C 2 (k) is determined by another boundary condition , C 2 (k) = 
co(k,0)/co2{k,0) , Then 


(o(k,z) = iku(0)h(k) 


tOn ( k 5 z ) 


Ck)- 


CO 


(X,2) = /”iku(0)h(k) 


,(k,0) 

033(k,z) 

G 

6o^(k,0) 


dk 


(25) 

(26) 


This integral is improperly defined for any value of k where ^^(kjO) = 0. 
If the integral is evaluated at these singularities by taking Cauchy's 
principal value, the primary contribution to the integral comes from the 
neighborhood of the singularities. These discrete values of k, if any 
exist, correspond to free waves or resonance waves of the system, and 
represent eigensol utions of equation (24) which dominate other waves in 
the system. 

For a given velocity and stability profile, the boundary conditions 
either do not uniquely determine the steady-state solution, or do not 
uniquely determine the amplitudes of the free waves. Mathematical 
uniqueness may be established, however, by some physical argument, such 
as requiring that all waves vanish far upstream of the barrier, or by 
considering a time-dependent system which asymptotically approaches the 
steady-state solution. 
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COMPARISON OF MODEL WITH ANALYTICAL SOLUTIONS 


To establish confidence in the consistency of the numerical model, 
the computations will be compared to known analytical solutions. In 
general, these analytical solutions exist only for the linearized steady- 
state equations with simple idealized meteorological profiles. Although 
the model also incorporates nonlinear and transient motions, it should 
qualitatively and quantitatively resemble the linear, steady-state solutions 
after a period of model time, assuming that the barrier height and the 
density and velocity profiles have been selected to preclude highly non- 
1 i near effects . 


Constant Velocity, Constant Stability Case 


This is basically the simplest case, since F(z) in equation (20) has 


the constant value F = \ where u = u , Usually, ^ s^« 

^0 Uq 

that the (:a' equation becomes: 


so 




(27) 


9x 


9z‘ 


,gS . 

where k z ( — ) is the stationary wavenumber. This wave equation specifies 

^ 2^u T 1/2 

a disturbance with wavelength A = -^ = - = 2ttu - P - - -- J (where T 

s N 0 qiYj -y] 

5 ^ ' d 


N 0 -qiY , ' 0 

is the reference temperature, y is the lapse rate, and Y^ is the dry adi- 
abatic lapse rate), which is independent of the barrier and is approached 
asymptotically at large distances by the actual solution since the boundary 
conditions have not yet been taken into account. Using the previously 
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stated boundary conditions, Lyra Cl 9^3) expressed the vertical velocity 
field for an arbitrary obstacle z = h(x) as an infinite series of Bessel 
functions of the first and second kind (J^ and Y^) : 


oi' (x ,z) = 2u^ /° 


dh(x) 

dx 


__8 

8z 


ci 


Vo(k,r) 


+ 1 

7T 


E 

v=0 




cos(2v+l )a 
2v+l 

(28) 


Jdx 


-1 X 

where a = tan {—) . Since the Bessel functions are eigensolutions of 
the system, an infinite number of free waves exist, due to the zeros in the 
J^. For many barrier shapes, including the rectangular shape in the 
numerical model , the free waves add up to form an infinite series of back- 
ward tilting lee waves with X as r ^ 

Figures 1 through 3 show the transient development of the streamlines 
for the Lyra problem when u^ = 25 m/sec, = 273K, y = 0, Ax = Az = 1000 m 
and At = 20 sec (referred to as case 1) at times 50 At, 75 At, and 100 At. 
The qualitative appearance of the waves agrees with the results of Lyra, 
except for the influence of the top boundary in the model. This boundary 
will not be as important in most other velocity profiles, since more 
energy will be trapped at lower levels, The theory predicts that X^ = 8.4km 
The most reliable wavelength measurements for comparison to the analytical 
solutions are taken in the area of a well -developed wave pattern and as far 
downstream as is feasible to avoid distortions caused by the assumption of 
no upstream perturbations, At 7 km elevation, figure 3 exhibits 8.5 km 
separation between the second and third crests, and 8 km separation between 
the third and fourth crests . 
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As the model approaches a steady state, 

^ / 1^; that is, the slopes of the isolines of streamfunction 

and density are everywhere equal, so that the isolines of these variables 
should coincide. Figures 4 and 5 show the general resemblance of the density 
field to the streamline field at time 50 At and 100 At, except in the vicinity 
of the barrier, where most transient development is still taking place. 

Assuming that p(z) s: p in the model, then co' ~ (— = w‘ = 

p 

and equation (20) becomes: 

2 , 2 

+ F(z)ip' = 0 (29) 

8x 9Z 

p 2 

In the Lyra case, F(z) = , and c = 0, so then e' = V i/j', which implies 

2 

that c' = -k^ that is, the isolines of vorticity should coincide with 
the isolines of streamfunction displacement. Figures 6 and 7 show the 
resemblance of the vorticity field to the troughs and crests of the 
streamline field, except in the vicinity of the barrier, where the 
assumption of linearity is not valid. It should also be noted that the 
vorticity field shows small scale perturbations which are computational 
in nature. These short wavelength perturbations occur mainly as a result 
of aliasing, that is, the inability of any numerical model to resolve 
disturbances with wavelengths less than two grid intervals. As is 
discussed in the appendix, the finite differencing scheme used retards the 
unstable growth of these perturbations, and we have not found it necessary 
to use filtering, smoothing, or damping operators in order to run physically 
meaningful computations for a sufficiently long period of time. Since the 
vorticity is the second derivative of the streamfunction, the streamfunction 
field should remain smoother than the vorticity field. 
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The local Richardson number, defined as Ri = - ^ in the 

P 32 / 32^ 

model corresponding to figures 3, 5, and 7 at time 100 At, is shown in 
figure 8. It should be noted that Richardson number tends to vary 
rapidly over several orders of magnitude in disturbed sections of the flow. 
Thus, log-|g(Ri) is actually plotted, and values of Ri < .16 or Ri > 10 are 
set to 0.16 or 10 respectively, in order to highlight areas where Ri < 0.25. 
Also, since the Richardson number is the quotient of first and second 
derivatives, the finite-difference analog for Ri is not dependable within 
one grid interval of the ground terrain, 


Linear Shear, Constant Stability Case 

S ^ 1 2 

For this case, F(z) = + where u = u^(l + cz) . Assuming 


u 


u 9z 


that >> ■§- equation (29) becomes: 

u u 


^ ^ + (4 - 4)r - 0 


(30) 


8x 


3z 


u 


CO 1 k X ^ / \ 

Assume that ip'(x,z) = / e i|j(k,z)dk, and let z^ = 1 + cz. Then, for a 

— oo 

single wave component: 

j2^ o Ri- ■- 

+ (-k/ + = 0 (31) 

dZi z^ 

k2^^ 

where k-. = ^ , and Ri = . The solution to this equation 

c ° u^c 

satisfying the upper boundary condition is a modified Bessel function of 


the third kind of imaginary order K. (k-jZ^), where y = (Ri^ - ^) 


1 . 1/2 
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Using the lower boundary condition, the solution can be expressed as (Wurtele, 


1953): 


, 1 1 \ ^l+cz^ 1/2,00 

^'(x,z) = (- 2 :jp) / 


ixk,^ 1^4,, 
h(k) 


[l <1 (Hcz )] 
^TkT) 


dk 


iy'"l 


(32) 


The free waves of the system correspond to discrete values of k^(or (ki)^) 
for which K^.^[(k-|)^] = 0, with wavelength = 27 r[(ki)^c^- j 

and no tilt with height, and exist only for p > 0 or Ri^ > The number 
of free waves with wavelengths in the mesoscale range increases with 
decreasing shear, approaching an infinite number in the Lyra case. 

When u^ = 1 0 m/sec, c = 2,7 x T^ = 250K, and y = 6,76K/km 

(referred to as case 2), then Ri^ = 16.0, p =3,97, and the theory predicts 
two free wave modes with wavelengths 13.7 km and 31.0 km. Figure 9 shows 
the streamfunction field at 3750 seconds. Only the first wave has developed, 
with two crests separated by 13 km. At 7500 seconds, figure 10 reveals the 
shorter wave dominating below 5 km, and a longer wave with two crests 
separated by 29 km prevailing at higher elevations, in accordance with the 
theory. 


Exponential Shear, Constant Stability Case 

T ~ cz gS -2cz , 2 1 2 , 1 3s 

In this case, u = u^e , so F(z) = e +sc - c - ^ s + ^ 

^0 _ 

Assuming that the total percentage variation of p is small compared to the 

2 9s 

total percentage variation of u, then the terms s , — , and sc are all much 

dZ 

2 

smaller than c , and equation (29) becomes: 

+ (^ - c^)i|j' = 0 (33) 

9x"^ 9z^ u 

0 
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Assume that ij;'(x, 2 ) = f° e^**^^Tp{k,z)dk, and let Z 2 = where 


Ri E — ^ , 

° u ,2 

° 2- - "T ^ 

^ + (1 = 0 

dZo 2 z^^ 


Then, for a single wave component: 


(34) 


The solution to this equation satisfying the upper boundary condition is 

k^ 1 /2 

a Bessel function of the first kind, where v = (—2 +1) ^ > 1. Using 

c'^ 

the lower boundary condition, the solution can be expressed as (Palm and 


Foldvik, 1960): 

1 ( Ri ^ 

V(x,z) .7(z)/" - 


dk 


(.35) 


The free waves, if any, result from discrete values of k^(or v^) for which 


J^(Ri^^^^) = 0, and have wavelength 


_ 2tt 


2tt 




with maximum 


amplitude development at z 


- — ln[(z„) /Ri ], where (z„) is the value 

V* ^JlV^ ^ll 

of 22 at which attains its maximum value between the nth and (n+l)th 

zeros of 3^(22^ • 

It can be seen graphically (Jahnke and Emde, 1945) that no waves exist 

1/2 1/2 
if Ri^ ' < 3,8, one wave exists if 3,8 < Ri^ ' < 7.0, and two waves exist 

if 7.0 < Rig^"^^ < 10,2. When u^ = 20 m/sec, c = l.SxlO'^^/m, and 

N = 1.2xl0~^/sec (referred to as case 3) , then =3.3, and the theory 

predicts no waves. Figure 11 shows no waves after 1200 seconds. When 

4 -2 

u^ = 10 m/sec, c = 2x10 /m, and N = 1 .2x10 /sec (referred to as case 4), 

1 /2 

then Ri^ ' = 6.0, and the theory predicts one wave mode with X = 12,8 km, 
and with maximum amplitude development at 2 = 2,24 km. At 4500 sec. Figure 12 
shows two crests separated by 13 km and maximum amplitude occurs at z =2 km. 
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Nonlinear Case with Constant pu'^ 

This case, developed by Long (.1955), solves the nonlinear equations with 
a nonlinear barrier in a fluid v/ith a rigid top and bottom. The incompres- 
sible, steady-state equations of motion (13a, b) may be rewritten as follows: 


2 2 
8 /U +w X 

P 3x ^ 2 ^ 

2 . 2 




3 /U+Wx^^ _ 3p 

( — 2 — ' ^P^ " " Fz " P^ 


3z 


(36a) 

(36b) 


where p = p(ip) and c = Eliminating p, equation (36) becomes: 

? I ^ "9Z)] = 0 (37) 


This is then integrated to yield: 

n2 , , 1 dp (Vii<) ^ , 1 dppU , f _xn 

V ^ -V- " ^ + p dt^-2 ^ S^^rx-Z)] 


P dTp 


(38) 


where u(ijj) and 2 ^( 4 ;) are the horizontal velocity and height associated with 


the streamline for constant ifj far upstream of the barrier. Noting that 

d\|j 

n 2 ^ 4 . riv 

2 Sz-' dz^ 

Now, if pu^ is constant, then 'is constant, and equation (39) 


5 and substituting 5 = z - z yields an equation in 6 : 

0 

„2„ ^ r(V5)^ ^ 36-, d ,, 2x g dp . 

V 5 + V o ■ + :tt Hn pu ) = -^2 ^ (39) 

'0 pu 0 


becomes the linear equation: 


2 2 
V 6 -I- a 5 =0 


(40) 


2 ci / 2 2 

where • ^'he constant pu criterion is satisfied approximately 

P 0 

by the model by setting u constant and keeping density gradients small. 
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Specifying the boundary conditions as 6(.x,0) = |-{1 + cos for 
-b < X < b, with 6(x,0) = 0 elsewhere, and 6(x,H) = 0 at the top boundary 
z = H, Long expressed the solution of (38) as: 


5(x>b, z)=-aZ [2Y„sin(5.pb}sin(ilrix) sin(nTTz)] 
n=l " 


- a 


f {In 

n=n^+l 


(x+b) 


]sin(n7Tz) } 


(41) 


2 2 2 2 riTT / 2 ^ 

where 5, = - n'^Tr'^, = — 2 ” (^n ~ and n-j is the largest integer 

2 ,2.n b 

for which 2,^ >0. 

When u^ = 30 m/sec, N = 1 ,2x1 0~^/sec , H = 10 km, a = 0.3H, and b = 0.4H 

(referred to as case 5), the theory predicts a single wave mode with wave- 

2 2 

length A-, = = 2TrH/(-— d = 25.4 km with maximum vertical velocity 

"'1 ' u 

0 

^\ax^ ~ ^^^0 'Y'l^iSin b = 21.2 m/sec. In order to obtain numerical 

computations correspondi ng to Long's solution, the model simulates the lower 

boundary condition by specifying the nonrigid flow boundary 
u air 

w(x,0) = sin(-^) for -b < x < b, and w(x,0) = 0 elsewhere. Values 

of a and b have been chosen to preclude highly nonlinear disturbances 
downstream. Figures 13 and 14 show the streamline and vertical velocity 
fields at time 3000 sec. The separation between the two crests is approxi- 
mately 25 km, and the maximum vertical velocity is approximately 18 m/sec. 

It should be noted that the vertical velocity is defined to have a greater 
value than this at the inflow, and that the measurements must be taken at 
least beyond the first crest. 
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Comparison of Nonlinear Effects 


In the cases above, the height of the barrier has been chosen to be 
sufficiently small so that the model will produce linear effects. However, 
features such as the reversed flow or rotors often observed in association 
with mountain lee waves arise from the nonlinearity of the barrier. This 
concept may be explained as follows, According to equation (22), the 
linearized surface boundary condition for the streamfunction may be written 

\|j'(x,0) = u^h(x) = u^hd{x) (42) 

where d(x) is a dimensionless profile function of order unity, and h is 
the maximum height of the barrier. In the constant velocity, constant 
stability case, for example, ijj'(x,z) is a function of k^h, so that if; 
may be expressed in the form: 


= tf; + ifj ' = -U Z «G( k X , k^z) 
0 s s 


(43) 


where G is of order unity, and G(k^x, 0) = d(x). The horizontal velocity 
is then: 


34 , , , 

- if ' % ■ li 


' — Him 


3G(k x,k z) 

S S ^ 


where 9G/3(k^z) is dimensionless and of order unity. The condition for 
reversed flow, or u < 0, is then approximately k^h > 1, Miles (1969) 
calculated the critical values of k h for flow over semi-elliptical barriers 
to be between 0.67 and 1.73, depending on the ellipticity of the barrier. 

Figure 15 displays the streamline field at 1500 seconds for the 


Lyra problem when k^h = 1.17 (referred to as case 6), showing that the 
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critical limit has just been exceeded, with reverse flow at some points 
in the field. This case is similar to case 1, except that Ax = Az = 750 m, 
and At = 15 sec. The corresponding density field is displayed in figure 16. 
A more highly nonlinear case with k^h = 1.95, T^ = 250K,.Ax = Az = 625 m. 

At = 10 sec (referred to as case 7) is portrayed in figures 17 through 19, 
at times 600 seconds, 800 seconds, and 1000 seconds, respectively . This 
sequence reveals the development of highly unstable configurations which 
break down realistically into rotorlike formations. It should be noted 
that the turbulence associated with the instability of the breaking wave 
is not simulated, The density field correspond! ng to streamline figure 19 
is shown in figure 20. The Richardson number fields corresponding to 
streamline figures 15 and 19 are shown in figures 21, and 22, respectively. 


COMPARISON OF MODEL WITH OBSERVATIONS 

Some of the most detailed observations of the atmospheric structure 
associated with mountain-induced waves have been taken over the eastern slope 
of the Rocky Mountains near Boulder and Denver. Boulder is located at the 
immediate base of the north-south range, and is susceptible to occasional 
downslope windstorms. On January 11, 1972, a particularly violent wind- 
storm with peak mean wind velocities of 30 m/sec, and gusts to 55 m/sec 
swept through the area. Lilly and Zipser (1972) derived cross sections of 
the potential temperature (figure 23) and horizontal wind velocity (figure 24) 
from aircraft observations taken during this event. The figures reveal a 
severe downslope windstorm and extensive mid-tropospheric clear air 
turbulence induced by a wave of large amplitude and wavelength. 
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Nonlinear numerical simulations of this case have been performed by Klemp 
and Lilly (1978), and by Peltier and Clark (1979), 

The wind and stability profiles are initialized for our model from the 
Grand Junction, Colorado sounding at OOOZ on 12 January 1972. Grand Junction 
is approximately 300 km upwind and at approximately the same elevation as 
Boulder. The lee slope of the Rocky Mountains in the vicinity of Boulder 
is reproduced as closely as is possible with a resolution of Ax = 2000 m 
and Az = 500 m. The upwind terrain is quite complex, but its mode! 
representation was not found to have an appreciable effect in the computation 
owing to partial upstream blocking. Since the potential temperature often 
varies on the surface of a high mountain, the density has been allowed to 
vary on the surface of the barrier for this computation. 

Figures 25 and 26 show the streamline field and the horizontal velocity 
field at 4250 seconds. The model reproduces many of the observed features 
of the mountain wave. The computed trough of the wave is located almost 
directly over Boulder, which is situated within one grid point of the lee 
slope. The computation shows the first crest of the wave to be 38 km 
downstream from the crest of the mountain, compared to an observed distance 
of 37 km. The wave shows no tilt with height up to the tropopause at 
approximately 11 km, then tilts back sharply into the stratosphere. 

In comparing the locations of the maximum and minimum wind velocities, 
it should be noted that figure 23 contains two sets of observations taken 
several hours apart. The winds in figure 24 are derived from the data taken 
during the time when the trough of the wave had moved over the lee slope of 
the mountain, presumably due to variation in the upstream wind or stability 
profiles, A study of numerous windstorms in the Boulder area by Brinkmann 
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(1974) reveals that the surface wind speed maximum is localized beneath the 
trough of the wave, and the output of the model is in agreement with this 
finding. 


SUMMARY AND CONCLUSIONS 

This report has described the development of a numerical model for the 
simulation of nonlinear, nonhydrostatic stratified flow over obstacles, This 
type of model is appropriate for the investigation of clear air turbulence 
associated with gravity waves resulting from flow over mountain ranges. To 
simplify the equations to be solved, the flow has been assumed to be two 
dimensional and incompressible, and the Boussinesq approximation has been 
made. These features have made possible a code which is sufficiently versa- 
tile and efficient to accommodate case studies using either idealized profiles 
or actual sounding data. 

Simplicity has also been retained in the boundary conditions. Distur- 
bances are generated by a rigid barrier, which is part of the lower boundary. 
The top boundary is also rigid, with periodic boundary conditions at the 
sides. Although these conditions require a sufficiently large computational 
field to produce physically useful results before contamination occurs, the 
computational speed of the program has always made this feasible. 

The consistency of the numerical model has been established by compar- 
ing the computations to known analytic solutions, and by comparison with 
mountain wave observations. The model reproduced qualitative and quant- 
itative features of the steady-state solutions, and also realistically 
simulated breaking wave patterns associated with highly nonlinear obstacles. 
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These tests provide confidence that the model may now be applied to 
observational data for further comparison. 
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APPENDIX: DESCRIPTION OF THE COMPUTER PROGRAM 


DESCRIPTION OF CALCULATIONS 


The System to be Solved 


The program solves the time-dependent system of equations 


9 ? _ 




gt 


If " 

2 

V = C 


(Ala) 

(Alb) 

(Ale) 


on a rectangular grid in the x-z plane, with rigid slip boundaries 
(tangential flow) at the top and bottom (where \p and p have constant, 
fixed values), and periodic boundary conditions at the sides such that 
the inflow at one side matches the outflow at the other. To facilitate 
the finite difference calculation of horizontal derivatives in the 
program, the second to last column is a duplicate of the first, and the 
last column is a duplicate of the second, Disturbances are generated in 
the flow by a rigid barrier of user-specified shape and size, which becomes 
part of the lower boundary. This system is pictured in the diagram below. 
Grid points in the x-z plane are indexed with the letters (i,j) beginning 
with (i,j) = (1,1). In this system, all variables will be assumed to have 
values only at the same discrete grid points (i,j). All finite difference 
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expressions are valid for all unique (i,j) in the grid unless otherwise 
stated. 



Several comments need to be made regarding the boundary conditions. 
First, since the streamfunction is held constant at all times along the 
lower surface and the barrier, the program is applying a nonlinear boundary 
condition in every problem. Thus, the potential exists for nonlinear 
features to form even when simulating a linear analytic solution. The dis- 
tinction between "linear" and "nonlinear" cases is determined by the degree 
of nonlinear behavior in the solution. 

Second, it should be noted that a periodic boundary condition in the 
x-direction allows disturbances which propagate to either lateral boundary 
to reenter through the opposite side, eventually contaminating the solution. 
The computational field must be given sufficient horizontal extent so that 
useful results are obtained before contamination occurs. This disadvantage 
is offset by the absence of reflection at the lateral boundaries, and by the 
flexibility of the model to simulate a wide variety of nonlinear flow prob- 
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lens without the necessity of devising a suitable set of open boundary con- 
ditions for each problem. 

Finally, we mention that the top boundary condition constitutes a 
rigid lid. Since this is highly reflective, the computational field must 
be given sufficient vertical extent so that useful results are obtained 
before significant reflection occurs. Due to the computational speed of 
the model, it has always been feasible to utilize a sufficiently large 
array for these purposes. 

Scheme for Solving the Equations 

The program uses an explicit, centered-time (leapfrog), centered-space 
finite differencing scheme with fixed boundary conditions on 4; and p to 
solve the system of equations (1). Assuming that ip, p, and ? are all known 
at time steps m -1 and m, then equations (la) and (lb) may be integrated to 
yield p and c at time step m +1 : 

- i |£]-2At + (A2a) 

PQ dX 

= J(p"’, c"’)-2At + p^"”* (A2b) 

where p"’ and are the values of p and c at time step m. Then is 

found from by relaxing . 

To begin the procedure, 4 ^, p, and ^ are initialized at all points at 
time step 0. Then, in order to start the three level time differencing 
scheme, it is necessary to obtain the values of the variables at time step 
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1 by another method. The model uses a Matsuno-type scheme, which involves 
taking a half time step forward, then using a centered time difference to 
proceed to time step 1, as follows: 


D°) - 1 IS- ] • ^ + ?° 


% 


^1/2 _ 1/ 0 ,Ov At , 0 

p = J(p , )• ^ + p 


relax for 


c'' = ^^-At + C° 

Pp 


)'' = J(p^/^, + p° 


relax for . 


(A3a) 

(A3b) 

(A3c) 

(A4a) 

(A4b) 

(A4c) 


This scheme results in less computational error, and is less destabilizing 
than taking a full forward time step. 

After hundreds of time steps, the leapfrog scheme may introduce a 
time-splitting instability into the solution of this nonlinear model. 

This instability may be suppressed by the occasional insertion of a time 
step made by a two level scheme (Mesinger and Arakawa, 1976), Thus, the 
program restarts the procedure every 25 time steps with the scheme de- 
scribed above. This is not the only method which may be used, but it is 
a standard numerical procedure (Arakawa and Lamb, 1977), 
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Initialization of ijj°, p°, and 

The user determines the unperturbed ij5°(z), ^(z), and E^(z) by speci- 
fying the initial horizontal velocity profile, u°(z), and either the initial 
stability profile, S°(z), or the initial temperature and pressure profiles, 
T°(z) and p°(z). In two dimensions, ijJ®{x,z) e - /[u°(x,z)dz - w°(x,z)dx]. 
When w°(x,z) = 0, and u°(x,z) is a function of z only, this becomes: 

^(z) = -/^u°(z')dz' + (A5) 

where is an arbitrary constant which has been set to zero by the program. 
Denoting the value of at grid point (i,j) as ij^., the finite difference 
form of equation (A5), using the trapezoidal rule, is: 


^,1 “ “ 


(A6a) 


tO _ 7-0 ( 

^i,0 ’^i.j-l 


1 >J-1 


— ■ )Az for all j=2,...,NJ (A6b) 


where ( j-1 )Az ] . 

^ 0 
For the compressible atmosphere, the stability is 5*^(2) = 

0 (z) 

The density profile in this incompressible model is defined to correspond 

1 ^—0 i 

to the stability of the compressible atmosphere by setting 
= S°(z). The expression for ^(z) is then: 


_^1_ 8^(z) 
—0 / \ 8z 

p {Z) 


^(z) = exp[-/^ s"'(z')dz'] 


rZ ^0, 


(A7) 
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3 

where is an arbitrary constant which has been set to 1 .25 kg/m . Using 
the trapezoidal rule, the finite difference form of equation (A7) is: 



(A8a) 


j exp[- + S°^j)] for all J=2 NJ (A8b) 


If the user specifies the temperature and pressure profiles instead of 

the stability, then using the definition 6° = , the density 

P 


profile is defined by setting 

“0 o Z ^ ^ 

for ^(z) becomes : 


R 3p 


c p z 
P 


The expression 


?^(z) = p. 


T°(z) 


P°(2) ^ 


(A9) 


where p^, T^, and p^ are arbitrary constants which have been set to 1,25 
3 5 2 

kg/m , 273K, and 10 kg/m-sec respecti vely , and k = 2/7. The finite 
difference form of equation (A9) is: 


^ / Pj ,j >2/7 

^i,j -j-o Pc 

i ,j 


(AlO) 


The vorticity is calculated initially from the streamfunction by the 
expression ^ Since ^ is a function of z only, the expression 

for ^ becomes: 


4? 


( 2 ) 


aV(z) 



(All) 
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The finite difference form of equation (All) is: 



(a1)' 



, 'tQ 

^ +i , j+1 



for j=2, , , . ,NJ-1 


(A12a) 



(A12b) 




,NJ-1 


(A12c) 


Influence of the Barrier 

2 

Since Fourier transform methods are used to relax V ip = ^ for the 
streamfunction at each time step, all interior grid points in the field, 
including those on or inside the barrier, must be considered to be part of 
the flow. Values of ip at these points are thus subject to change as a re- 
sult of the transformations. To preserve the desired boundary condition on 
the barrier, the effect of the vorticity generated by each separate point 
on the barrier is superposed with the effect of the vorticity generated by 
all the internal points in the grid. Then, the streamfunction at each point 
on the barrier is expressed as a linear combination of the relaxation 
solutions associated with these vorticities (Roache, 1972): 


ip^^ = ip^ + E a.ip^ for all £ = 1, NPB (A13) 

0 k=l K K 

2 2 
where \p^ is the solution of V \p^ = ip^ is the solution of V 

due to a unit vorticity at barrier point k, the superscript £ represents 

the values of and at barrier point £, and NPB is the number of 
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points on the barrier. The are determined at each time step from the 
linear system (A13), Then, at each grid point (i,j) in the system. 


NPB 




(A14) 


Although the superposed solution ^ results from additional vorticity on the 

2 

barrier, the solution is a valid one, since it satisfies V ij; = ?; at all 
internal points in the flow, and satisfies all boundary conditions, includ- 
ing the barrier. The additional vorticity on the barrier introduces no 
perturbations in the vorticity of the flow at any time. 

As an example in computing the ctj^'s, consider the case with three points 
on the barrier. The system of equations (A13) then becomes; 


4 ''* = 'pQ “l^l “3^3 


2 2 2 2 2 
^ + ct2^2 “3^3 


.3 


3 3 3 3 

ii ^ + “3^3 


Applying Cramer's rule. 


“1 D 








^0^ 


.^2) 


- 4^^) 
^0 ’ 


,1 

^^2 


.2 

^2 


|3 

^2 


,1 

^3 


,2 

^3 




(A15) 


(A16) 
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where D is the determinant of the matrix 


’'•’1 

,2 

^1 

r3 


,1 

^2 

,2 

,3 

^2 


,1 

^3 


V'. 


and B 


n D 


1 

,2 

,2 

*3 

a 1 

,1 

4>2 

,1 

'*'3 

o 1 

,1 

4^2 

,1 

*3 

D" 

,3 

,3 

’ ^12 D 

,3 

,3 

’ ^13 D 

,2 

,2 


\l>2 



4^2 

*3 


4^2 



(A17) 


with similar expressions for a,^ and . For the general case with NPB 
points on the barrier. 


a,^ = Z ^ 


(A18) 


where the need to be calculated only once. A Gaussian elimination 
scheme is used to calculate the determinants for the 3|^£'s. 


Specifying the Initial Conditions 

At time step 0, the barrier is suddenly introduced into the flow by 

defining the bottom topography to be a line of constant \p and p. To reduce 

the physical shock of introducing the barrier, a solution due to the barrier 

may be added to without perturbing the vorticity in the flow. This is 

expressed as + 4^barrier’ streamfunction at 

2 

time step 0, and V ^ 0- Since the isolines of i|j and p coincide 

everywhere as the model approaches a steady state, a similar perturbation 


is added to ^ so that the isolines of and p° approximately coincide, as 
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shown below. 


^ 4 ^ -5 > P 


P = P3 
- 4^2 



The finite difference expressions for ip° , p°, and are then: 


^ ^ ^ NPB 

4^? . = 4^ . + z {(4 jJ^ ,• [4^„ - (/)°]> 

> 5J 1 Sj |^_1 ^ I SJ ^ 


(A19a) 


0 _ — 0 / 1 0 tP \ ^ — P 

pi.j - pi,j'-i ^ 


~P \ Jf~rO 7-0 \ 

(Al9b) 


^P _ 


(A19c) 


where, in equation (19b), j' >2 is the lowest row number for which {ijn .,| 

5 J 

It should be noted that the potential flow scheme described above does 
not have a beneficial effect for every possible initial velocity profile 
u°(z). Specifically, if the velocity changes direction at upper levels, or 
is generally decreasing with height, then the code should be modified to 
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set \fj^ and p° to and (except on the barrier, where ip^ = and p° e p^). 

Calculation of J(p,ijj), and 

aX 

The long-term computational stability of the model depends on the finite 
difference form of the equations to be integrated. Arakawa (1966) devised 
a method to retard nonlinear computational instability in the equation 
by conserving mean vortidty, mean kinetic energy, and mean 
square vorticity. This scheme is used by the program to calculate J(i;»^|^) on 
the boundaries, and and J(p,i|j) at interior grid points. The finite 

difference form, J^.j(^,\j)), depends on the location of the grid point, which 
may be one of ten types, as shown below. 


7 (left outside 
corner) 


3 (left side 
of barrier) 


0 (interior) 


2 (top of 
\^^rri er) 






w 

L 


1 (top) 


t 


8 (right outside 
corner ) 


4 (right side 
of barrier) 




2 (bottom) 5 (left inside 9 (inside ’6 (right inside 
corner) barrier) corner) 


After applying boundary conditions, and using the fact that Is constant 
on the upper boundary and e 0 on the lower boundary, the finite difference 
expressions for J(c,\p) are as follows: 
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-1 


for type 0, ^+1 ,j-l " " ^+1 


^^"^i-l,j-l ' ^i,j-l ’*'^i-1,j+l ^ ,j+l ^^1-1 ,j 




^ ,j+l ^^i+1 , j+1 ^ ^ ’^i-l ,j^^i-l ,j-l 




-1 


for type 1 , J. .(^,Ui) = v~ ■ .- [(;ij. . , + iJj. ,n . n - 2\b. .)r. ,-, 
j-K . 6 AxAz ^1+1 ,j-1 1+1,; 




^■'^i+1,0-1 ■*■ '^i-l,j-l^^i,j-l 


’*’ ^“^i,j-l '*’ ,j^^^i-l ,j-l ' ^i+l,j-l^^ 


for type 2, J+] ^ j 


^"^i-l,j+l ’^1 ,j+l ^'°i-l ,j 


^'^i+l,j+l " ^1-1 ,j+l ,j+l ,j+l ^^1-1 ,j+l ■ ^i+l,j+l^J 
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,J+l^^i+l ,j 


(A20a) 


(A20b) 


(A20c) 




II 



+ J+l)] 

n -1 /_.\ -1 r/. 


i j 


"^i+1 J^^i+1 J+1 " ^i+1 


ij 


-1 ,j-l '^i-l ,j^^i ,j“'' 

(A20d) 

" ’^i+l J+l^^i+l,j 
•1 ,j-l " '^i+1 ,j^^i ,j-l 

(A20e) 

1-1 ,j ” ^i ,j+l 

(A20f) 

i,j+l ■ ^i+l,j^^ 

(A20g) 


for type?, J..(c,ip) = ^ , j 

^ ^"^i-1 ,j-1 ^^i-l,j+l ^ , j+1 ^^i-1 , j ^ 

'^i-l,j “ "^i-1 J+l^^i,j+l "^i-1 

for type 8, J.j{5,4,) = ff^+l ,j.l - l-i j +1 ‘ Vij+i )C ,+1 ,j 

^ ,j+l ^ ^1 J+1 ^^i-1 ,j ^ ,j ^ '^i+l ,j+l " ^i-1 ,j+1 ,j+1 

^“'^i+l,j-l " '^i+l ■" ’^i+1 ,j+l " ^i+l,j-l^ 

"^i,j+l ^^i-1 ,j+l " ^i+l,j+l^^ 


for type 9, J^.^.(c,^) = 0 


(A20j) 


For point type 0, J..(p,Tj;) has a similar expression to equation (A20a) 

* u 
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For all other types, = 0, which is equivalent to fixing the value 

of p on the boundaries. This method appears to be the most suitable for 
modeling physically significant transient motions in the system. When the 
Arakawa form of the Jacobian was used for J(p,ijj) on the boundaries, the 
system remained stable, but often approached a steady state at an unaccepta 
bly slow rate, since p was free to vary on the boundaries while was fixed 
The finite difference expressions for are: 


= ^1+1 .j ~ ^i-1 J 

^8x'i,j 2 Ax 


for point types 0,7,8 


(A21a) 


= 0 


= ~ ^i-lJ 

Ax 



for types 1 ,2, 5, 6, 9 


for type 3 


for type 4, 


(A21b) 

(.A21c) 

(A21d) 


Relaxation of V‘^ip = ? for ip 


The finite difference form of the equation V ip = <; is 


■1 ,J 


(Ax)^ 


(Az)^ 


(A22) 


The exact, noniterative solution of equation (A22), with periodic boundary 
conditions in x, may be expressed as the sum of discrete Fourier components 
c^(z) in X at each level of z: 
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(A23) 


Nx-1 

<.(X,Z) = i X C (7)6*"'^’' 
X n=0 


where k = ; x = pAx for p=0,. . N^-1 ; z = qAz for q=l , . . N^; 

X 

= NI - 2, the number of unique grid points in the x-direction; and 
= NJ - 2, the number of interior grid points in the z-direction. 


The c^(z) are calculated from the boundary values of ip, and from the Fourier 
components d^(z) of c (where c(x,z) = ^ d^(z)e^'^'^^) , through a system 

of finite difference equations which result from marching in the z-direction 
for each component n. 

The relationship between the Fourier components c^ and d^ is derived as 
fol 1 ows ; 


?(pAx,qAz) = 


Nx~l 

1 Z d (qAz)e 
n=0 ^ 


inkpAx 




(pAx ,qAz )^ 9 ip(pAx,qAz) 
3(pAx) 9{qAz)^ 


= i!r ( 


Nx-1 


''^x (Ax)^ 9p^ (Az)^ 9q^ n=0 


) E c (qAz)e^'"^P^^ 


(A24) 


= i e'"\c (qAz) ^ c (qAz)] . 

^x n=0 (Ax)^ dp^ (Az)^ dq“^ 


Using a centered finite difference approximation for the second derivative. 


1 ^ JnkpAx ^ 

(Ax)^ 


(Ax )^ dp^ 


(A25) 


= einkpAx _2 [cos(nkAx)-! ] 


T .2 C [(q+l)Az] + c [(q-l)Az] - 2c (qAz) 

and — L_^ c^(qAz) = 


(Az)" dq‘ 


(Az)‘ 


(A26) 
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Equation (A24) then becomes; 


Nx-1 


Nx-1 


^ = Z {c_ [cos(|l^)-l]e^‘"^P^>^ 


n=0 


n,q 


n=0 


n,q 


(Ax)‘ 


. ^n,q+l ^ ^n,q-1 ' ^^n,q ^inkpAx 
+ P e 

(A2)2 


} 


(A27) 


where ^ = c^{qA 2 ), Since each Fourier component responds independently 
of the others, the relationship under the summation holds separately for 
each n: 


] = c 

n,q n,q 


- 2c 


[cos(-P) - 1] + (A28) 

(Ax)^ ^ (Az)^ 


which becomes: 


-c 1 + [2 + (k')^]c ^ - c ,, = (Az)^d 

n,q-l \ i- J n,q n,q+l ' ' n,q 


(A29) 


where {k')^ = 2 [1 - cos(|S!-)], 

(Ax)^ "x 


For each n=0,. . , , - 1 , equation (A29) comprises a set of 

linear equations in N unknowns c , for q = 1,, , .. N , and is equivalent 

z n 5 q z 


to the matrix equation Ac^ = b: 
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-1 


-1 2+k' 


-1 2+k' 


-1 2+k' V \ c / 

/ N^.N^X "‘V N^.l 


“„,1 


-(Az)‘d 


n,N_-l 


“^n.N+l - >^n,l 


V Nz.l 


(A30) 


The d terms are found by taking the Fourier series of r at each interior 

j H ^ 

level z = t|AZ, q = 1,. , N , and the c and c ^ terms are found by 

^ rijD 

taking the Fourier series of at the upper and lower boundaries, as follows 


_ = 2 ^(pAx,qAz)e 

p=0 


-inkpAx 


(A31a) 


- = E i|j(PAX,0)e 
p=0 


-inkpAx 


N^-1 

^n N +1 = ^ i|^[pAx,(N +l)Az]e"'"*"P^'' 
’ Z ' p=0 ^ 


(A31b) 


(A31c) 


From this, the c terms for each n = 0,. . N - 1, for all q=l , , , .,N. 

* * * H X i 


are found by inverting the tridiagonal matrix A: 


= A~^b, for each n=0,...,N^-l . 


(A32) 


Then, ijj is obtained by taking the inverse Fourier transform of the c 


n,q 


N -1 


1 


i|;(pAx,qAz) = ^ E c e 
X n=0 


inkpAx 


(A33) 


The scheme for calculating the c^ ^ for each n is the following: 


a = 2{1+ [i_cos(f^)]> 

(Ax)^ \ 


u (1) = aj, f(l) = b(l) 


then u„ = - — 'r—T\ 

I Up(q-l) 


Up(q) = 


f(q) = b(q) - u^- f(q-l ) 


V for q=2, , . . ,N. 


(A34) 


then c 


f(N^) 
n,Nz = yV 


and c. 


. S.q+1 


n,q 




for q=Nz"l , . . . ,1 
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Calculation of u, Ri , and w 


These quantities are calculated from ii, p, and c at any time step 
specified by the user. Since the values are not used for any subsequent 
time step, the finite difference equations for u are smoothed by using cubic 
spline function interpolation to calculate continuous first and second 
derivatives of ij; and p in z. 

For each column of x, the procedure is to match the first and second 

derivatives of subsequent pairs of NJ-1 cubic polynomials q.(z) (j = 2,.,,, 

3 

NJ), at NJ-2 internal grid points, z . ( j=2 , . , . ,NJ-1 ) , given boundary condi- 

J 

tions at and Let these polynomials have the form: 


qj(z) = tVj + (l-t)Vj_^ + tAz(l-t)[(kj_^ - dj)(l-t) - (kj-dj)t], (A35) 


z-z 


where t = 




Az 




-VLh 


Az 


, V - = the value of or p at point j, and 


dq-(z J 


3 3 


(Dahlquist and Bjorck, 1974). Then the expressions for the 


^3 dz 

first and second derivatives in z become: 


dz 


+ (3t^-4t+l)(kj_^-dj) + (3t^-2t)(kj.-dj) 


dq . (z) v .-V • 1 
^3 ^ = 3 3 -1 

d 

and — ^ [(6t-4)(k^._^-d^.) + (6t-2)(kj-dj)] . 


(A3 6) 
(A37) 


dz' 


These polynomials satisfy the relationships: 




dz 


dz 


= k. 
3 


for j=2 NJ-1 


(A38) 
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and 


d^q.(z.) d^q ,(z.) , 

■ ^2 ^ j=2,.,..NJ-l (A39) 


provided that • 


(A40) 


The boundary conditions are: 


2 2 

^qp(^i) SV’I •! 

^2^ = " A? ^^2 “ ^*^2^ 2k^+k2=3d2~ ^ 

dz 8z 


a 'AZ 




(A41a) 


dz‘ 


3Z 




+ 2k., , = 3d., , + 


b'Az 


'NJ-1 NJ '""'NJ 2 


(A41b) 


Equations (A41a), (A40), and (A41b) comprise a set of NJ linear equations in 

NJ unknowns k., j=l,,..,NJ, and are equivalent to the matrix equation Ac =b: 
J ^ 



1 0 
4 1 

1 4 ’ • 


'. 4 

1 

0 




(A42) 


NJ,1 


3v . 

The iT“ = k. are found by inverting the tridiagonal matrix A according to 

oZ j 

the following scheme: 
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II? 


Up(l) =a^(l) 


f(l) = b(1) 


■1 


" Up(j-1 ) 

u„(j) = a ,(j) + u, 


f(j) = b(j) + u^.f(j-l) 


A for j = 2 , . . . , NJ 


(A43) 


then k,M , = 

NJ Up{NJ) 

f(j) - kj+1 for j = NJ - 


where aj(t)=2, a^(j)==4 for j=2 NJ-1 , and a^(NJ)=2. 

Then, by equation (39) the second derivative of v. becomes: 

J 


9^v 


9z‘ 


- = -i C-'* ^ ^ NJ-' 


(A44) 


The finite difference expressions for u = 


are as foil ows : 


-§ and Ri 


/ 

p 9z ' 


9z‘ 


,3 




9z 




-g 


9p 


9^1^, . 


i,j 9z 


^ 

z ' Sz 


(A45) 

(A46) 


The first derivative of is obtained for each column (i=l,,..,NI) by 

i~^i i-1 - -n 

applying equation (A42) with d. .= — . ■■■ — , a'= ? . , and 

1 5 J A Z 1 5 ^ 

b' ^ ^ Then the second derivative of ijj may be expressed as: 
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.J 

3z^ 


- J r_4 1 >J + g/ 1 a ‘ »J \ _ 

Az '- 9z ^ Az ^ 


3^. -^1 
3z -* 


(A47) 


The first derivative of p is obtained for each column (i=1,...,NI) 
by applying equation (A42) with d . = 

1 5 J AZ 


Pi i “ Pf i_i 

= _LiJ , a' - 0, and b' 


= 0 , 


The finite difference expression for w = ^ is: 

d X 


. _ '^i+l J " ^i-1 ,j 

m 


(A48) 


This equation may be solved by expressing Tp(x,z) and w(x,z) as a sum of 
discrete Fourier components c^(z) and dj^(z) in x at each level of z; 


N^-1 

iHx,z)=j1e C(2)e'"'"' 

n=0 

w(x,z)=^Z d (z)e''"'''^ 
X n=0 " 


(A49) 

(A50) 


2-it 

where k = , x = pAx for p=0,,,,,N -1, and N =NI-2, The rela- 

A 

tionship between the Fourier components c^ and d^^ is derived as follows 


X n=0 


N -1 

1 ^ 

— — — E c (z)e 

N^AX 3p n 


N -1 

^ inkpAx _ _L 1 fL ^inkpAx 


^ n=0 '’P " 


N -1 

X 


, A ink(p + 1 )Ax -ink(p - 1 )Ax 

“ r -i 5in(nkAx)e'’"''P*>' 

X n=0 


(A51) 
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Since each Fourier component responds independently of the others, the 
terms under the summation are equal for each n: 


d^(z) = c^(z) ^ sin(nkAx) (A52) 

The c^(z) terms are fo*und by taking the Fourier series of ip at each 
level of z: 

c (z) = E i|;(pAx,z)e"'"^'^P^^ (A53) 

^ n=0 

Then equation (A52) yields d^(z), and w is obtained by taking the inverse 
Fourier transform of the d (z): 

w(pAx,z)=^E d (z)e^"'^P^^ (A54) 

n=0 
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1 1 ■■■■ 


DESCRIPTION OF CARD INPUT DATA 
Space and Time Data . 

Space and Time Card : N I , N J , NT , NPB , I BT , I SAVE , DELTAX, DELTAZ , DELTAT 

(format: 614,1 P3D1 0.2) 

NI Number of grid points in the x-direct'ion (i.e., number of columns). 

Due to restrictions imposed by the fast Fourier transform 
routine, NI must equal 2 ^ + 2, where n is a positive integer. 

NJ Number of grid points in the z-direction (i.e., number of rows). 

NT Number of time steps. 

NPB Number of grid points on the surface of the barrier, excluding 

points on the lowest row. Refer to explanation of barrier 
data below. 

IBT Beginning time step. If IBT = 0, then read wind velocity and 

stability data. If IBT.NE.O, then the data saved at time 
step IBT from a previous run is to be input from a tape or 
permanent file. 

ISAVE Data retention indicator. If ISAVE = 1, then data at the last 

time step is to be saved on a tape or permanent file. Other- 
wise, set ISAVE = 0. 

DELTAX Grid spacing in the x-direction (in meters). 

DELTAZ Grid spacing in the z-direction (in meters). 

DELTAT Time step interval (in seconds). 
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Barrier Data 


Barrier Card(s) : (IB(I) , JB(l) ,1=1 ,NPB) 

(format: 2014) 

IB(I) Column number of the Ith point on the barrier. 

JB(I) Row number of the Ith point on the barrier. 

As shown in the example below, IB(I) and JB(I) are to be specified in a 
continuous fashion along the surface of the barrier, beginning with the 
leftmost grid point in the second row, and ending with the rightmost grid 
point in the second row. The surface of the barrier must connect adjacent 
grid points only horizontally or vertically, never diagonally. No barrier 
grid point may ever be specified in the first row. The position and shape 
of the barrier may be arbitrary, but only. one barrier is permitted, and its 
width may not exceed (NI-4) horizontal grid intervals. On the barrier 
pictured below, NPB = 16, with IB (1) =5 JB(1 ) = 2, IB(2) = 6, JB(2) = 2, 
etc. The width of the barrier is 7 horizontal grid intervals. 
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Graphical Output Data 


The output of the program consists of shaded, line printer graphical 
displays of the streamfunction (tP), density (p), vorticity (?), horizontal 
velocity (u), vertical velocity (w), or Richardson number (Ri) at any time 
step specified by the user. The shading consists of alternating clear areas 
and printed areas. To emphasize aspects of the flow, the user may output 
any segment of the entire grid, vary the vertical scale of the printout, or 
vary the resolution of the shading. 

Grid Plot Card : IS, JS, IN, JN,NJSMl 

(format: 514) 

IS First grid point in the x-direction to be plotted. 

JS First grid point in the z-direction to be plotted. 

IN Number of grid points in the x-direction to be plotted. 

JN Number of grid points in the z-direction to be plotted. 

NJSMl Number of lines of print between two vertical grid points. 

Shading Level Card : NPSILV,NRHOLV,NZTALV,NULV,NWLV,NRILV 

(format: 614) 

NPSILV Maximum number of levels of shading on ip graphs. 

NRHOLV Maximum number of levels of shading on p graphs. 

NZTALV Maximum number of levels of shading on z graphs. 

NULV Maximum number of levels of shading on u graphs. 

NWLV Maximum number of levels of shading on w graphs. 

NRILV Maximum number of levels of shading on Ri graphs. 
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Plot Card : 


IPSIGR 
p Plot Card : 

IRHOGR 

g Plot Card : 
IZTAGR 

u Plot Card : 
lUGR 

w Plot Card: 
IWGR 

Ri Plot Card : 
IRIGR 


IPSIGR(I) (I<20) 

(format: 2014) 

Sequence of time steps at which rjj is to be printed. 

IRHOGR(I) (I<20) 

(format: 2014) 

Sequence of time steps at which p is to be printed. 

IZTAGR(I) (I<20) 

(format: 2014) 

Sequence of time steps at which is to be printed. 

lUGR(l) (I<20) 

(format: 2014) 

Sequence of time steps at which u is to be printed. 

IWGR(I) (I<20) 

(format: 2014) 

Sequence of time steps at which w is to be printed. 

IRIGR(I) (I<20) 

(format: 2014) 

Sequence of time steps at which log.|Q(Ri) is to be printed. 
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Note that all variables need not be output at the same time steps, but that 
all time steps must be in order for any one variable. If a particular 
variable is not to be printed at all during a job, then a "-1" must be 
punched in columns 3 and 4 of the appropriate card. Time step 0 is a 
legitimate time step at which any field may be printed. 

If IBT.NE.O, there are no more cards to be read beyond this point . 

Wind Velocity Data 

This data determines the initial values of ijj^(z) and ^(z). 

Wind Card : ICASE 

(format: 14) 

ICASE Wind velocity profile indicator. 

If ICASE = 1, this is a sounding data case, and temperature and pressure 
data are to be read in addition to wind data. 

Sounding Cards : U( J) ,T( J) ,P( J) for J = 1,...,NJ 

(format: 1P3D10.2) 

U(J) Horizontal wind at row J (in knots) 

T(J) Temperature at row J (in °C). 

P(J) Pressure at row J (in mb). 

If ICASE = 1, there are no more cards to be read beyond this point. 


If ICASE = 2, this is a constant velocity profile case. 



Velocity Card : U1 

(format: IPDIO.2) 

U1 Horizontal velocity at all levels (in m/sec). 

Then, U(J) = U1 for J = 1,... ,NJ. 

If ICASE = 3, this is a constant shear profile case. 

Velocity and Shear Card : Ul ,C 

(format: 1P2D10.2) 

Ul Horizontal velocity at row 1 (in m/sec). 

C Vertical shear (in sec~^). 

Then, U(J) = Ul + C * DELTAZ * (j - L) for J = 1,...,NJ. 

If ICASE = 4, this is an exponential profile case. 

Velocity and Shear Card : U1,U2,C 

(format; 1P3D10.2) 

Ul Constant velocity to be added to the profile at all 

levels (in m/sec) . 

U2 Base horizontal velocity (in m/sec). 

C Vertical, shear divided by U2 (in m~^ ) 

Then, U(J) = Ul + U2 * EXP(C*DELTAZ*( J-L) ) for J = 1,...,NJ. 

If ICASE = 5, this is a hyperbolic tangential profile case. 

Ve locity Card : Ul , U2 , L I , MDZ 
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(format; 1 P2D1 0.2,214) 


U1 Horizontal velocity to be added to the profile at 

all levels (in m/sec). 

(J2 Base horizontal velocity (in m/sec). 

LI Row number at which the profile has its inflection 

point. 

MDZ Number of rows away from the inflection point at 

which in is deflected by the amount U2. The sign 
of MDZ determines the sign of (U{NJ)-U(1 ) ) . 

Then, U(J) = U1 + U2 * DTANH((J-LI )/MDZ) for J = 1,...,NJ. 

If ICASE = 6, this is a case, other than a sounding data case, for which 
the velocity is to be specified at each row. 

Velocity Card(s) : U(J) for J = 1,...,NJ 

(format: 1P8D10.2) 

U(j) Horizontal velocity at row J (in m/sec). 

Stability Data 

This data determines the initial values of ^(z). 

Stability Card : JCASE 

(format: 14) 

JCASE = Stability profile indicator. 
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IF JCASE = 1, this is a constant lapse rate case. 


Lapse Rate Card : GAMMA, TO 

(format: 1P2D10.2) 

GAMMA Lapse rate (in °C/m). 

TO Reference temperature (in °C) 

Then, S(J) = (DALR - GAMMA) / TO for J = 

where DALR = dry adiabatic lapse rate. 

If JCASE = 2, this is a constant Brunt-Vaisala frequency case. 

Frequency Card : BV 

(format: IPDIO.2) 

BV Brunt-Vaisala frequency (in sec”^ ) . 

Then, S(J) = BV ** 2 / G for J = 1 , . . . ,NJ ; where G 
acceleration of gravity. 

If JCASE = 3, this is a constant Richardson number case. 

Ri Card : RI 

(format: IPDIO.2) 

RI Richardson number. 

Then, S(l) = R * (U(l) - U(0)) ** 2, where 
R = RI / (G * DELTAZ ** 2), 

S(J) = R * (U(J+1) - U(J-D) **2/4 



for J = 2, . . . ,NJ-1 , and 

S(NJ) = R * (U(NJ) - U(NJ-l)) **2. 
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Summary of input card data 












K. read ip, p, z, 

3 . 

ITYPE, a^, and 
or disk 


















LIST OF VARIABLES IN THE PROGRAM 


For reference, the variables are briefly described in alphabetical 
order, with array dimensions, and symbolic references, if any. 


A(NI,NJ) 

AD(NI) 

AJZP 

AK(NI) 

ARGI 

AS(NJ) 

A1 through A9 

BETA(NPB,NPB) 

BLANK4 

BV 

B1 through B8 
C 

CDZ 

CFAC 

CKTMPS 

CRI 


Variables in the Main Routine 

Temporary array for storage of fields to be plotted. 

2 

Diagonal elements of the matrix used to relax V ip = ^ 
for X. (a^) 

Jacobian for <; and \p, J(4,^) 

Coefficients used to calculate vertical velocity from ip. 
Argument used for sinusoidal functions in AD and AK. 
Diagonal elements of the spline matrix, (a^) 

Temporary storage for values of \p at selected grid points. 

Coefficients used to calculate a. (3) 

4 blank print symbols. 

Brunt-Vaisala frequency. (N) 

Temporary terms for the calculation of Jacobian terms. 

Shear in horizontal velocity, (c) 

C divided by DELTAZ. 

Cofactor sign divided by the determinant of the a matrix. 
Conversion factor from knots to meters/sec. 

Constant Richardson number. 


DALR 

DELTAT 

DELTAX 

DELTAZ 

DMGDTO 

DRHO 


Dry adiabatic lapse rate, (y^) 

Time step interval, (At) 

Grid interval in the horizontal direction. (Ax) 
Grid interval in the vertical direction. (Az) 
(DALR - GAMMA) / STMP, 

Density differential in the horizontal. 
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DXM2 

DXSQ 

DZD2 

DZD6 

DZM2 

DZSQ 

DZ2DX2 


DELTAX multiplied by 2, 
DELTAX squared, 

DELTAZ divided by 2. 
DELTAZ multiplied by 6, 
DELTAZ multiplied by 2. 
DELTAZ squared, 

DZSQ divided by DXSQ, 


FAC 
FACM2 
FACM4 
FACM43 
FTI(NI ,NJ) 
FTR(NI,NJ) 
FI (NI) 
F2(NI ) 


-1 / (12 * DELTAX * DELTAZ) 

FAC multiplied by 2, 

FAC multiplied by 4, 

FAC multiplied by 4/3. 

Imaginary part of the Fourier transform of Tpf (Im(C )) 

n 5 c] 

Real part of the Fourier transform of ip . (Re(c )) 

n 5 C) 

utility array used by subrouting CALCW. 

Utility array used by subroutine CALCW, 


G Acceleration of gravity, (g) 

GAMMA Lapse rate, (y) 


I 

IB(NPBP2) 
IBT 
IB25 
I CASE 
IDUM(8) 

IE 

lEB 

lET 

lEW 

IEXP2 

IL 


An index. 

Column numbers of points on the barrier. 

Number of the first time step. 

Indicator of an even multiple of 25 time steps, 

Wind velocity profile indicator. 

Filler array for tape or disk read/write. 

Number of the last column to be plotted. 

Rightmost column number on the barrier. 

Number of the last time step. 

Number of the last column at which w is calculated. 

Power of two of the number of unique horizontal grid 
points, N^. 

Number of print symbols to be plotted in a line of graphic 
output . 
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IN 

INMl 

IPSIGR(20) 

IRH0GR(20) 

IRIGR(20) 

IS 

ISAVE 

ISB 

ISW 

ITYPE(NI,NJ) 

IUGR(20) 

IWB 

IWGR(20) 

IZTAGR(20) 

J 

JB(NPBP2) 

JCASE 

JE 

JMAX 

JMAXMl 

JN 

JS 

K 

KM2 

L 

LI 

LINE14(33) 
LINE24(33) 
LINE3(132,NJSM1 ) 
LINE34(33,NJSM1 ) 


Number of columns to be plotted. 

IN minus 1 . 

Time step numbers at which is plotted. 

Time step numbers at which p is plotted. 

Time step numbers at which log-|Q(Ri) is plotted. 

Number of the first column to be plotted. 

Tape save parameter. 

Leftmost column number on the barrier. 

Number of the first column at which w is calculated. 

Grid point type. 

Time step numbers at which u is plotted. 

Horizontal grid point width of the barrier. 

Time step numbers at which w is plotted. 

Time step numbers at which c is plotted. 

An index 

Row numbers of points on the barrier. 

Stability profile indicator. 

Number of the highest row to be plotted. 

Highest row number on the barrier. 

JMAX mi nus 1 . 

Number of rows to be plotted. 

Number of the lowest row to be plotted. 

An index 
K minus 2. 

An index. 

Row number at which the hyperbolic tangential profile 
has its inflection point. 

Part of annotation on graph of field. 

Line of print symbols associated with a level of the grid. 
Logical symbol equivalent of LINE34. 

Lines of print symbols between two levels of the grid. 
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LINE4804) 

LINE58(14) 

LINE68(14) 

M 

MDZ 

MT 

N 

NAME(2,6) 

NEWRHO(NI,NJ) 

NEWZTA(NI,NJ) 

NI 

NIMl 

NIM2 

NIPA 

NIS 

NISMl 

NJ 

NJMl 

NJM2 

NJSMl 

NPB 

NPBMl 

NPBPl 

NPBP2 

NPSIGR 

NPSILV 

NRHOGR 


Part of annotation on shading value scale. 

Part of annotation on shading value scale. 

Part of annotation on shading value scale. 

An index. 

Number of rows away from the inflection point at which 
Ul is deflected by the amount U2. 

Time step number. 

An index. 

Titles of quantities on graphical output. 

Temporary array for newly calculated density. 

Temporary array for newly calculated vorticity. 

Number of grid points in the horizontal direction, or 
number of col umns . 

NI minus 1 . 

Number of unique columns = NI minus 2. 

Number of unique columns, plus the horizontal grid point 
width of the barrier. 

Number of print symbols from one horizontal grid point 
to the next. 

NIS minus 1 . 

Number of grid points in the vertical direction, or 
number of rows , 

NJ mi nus 1 . 

NJ mi nus 2 . (N^) 

Number of print lines between two vertical grid points. 
Number of grid points on the surface of the barrier 
NPB minus 1 . 

NPB pi us 1 . 

NPB plus 2. 

Next time step number at which ij; is to be plotted. 
Maximum number of levels of shading for graphs. 

Next time step number at which p is plotted. 
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NRHOLV 

NRILV 

NT 

NUGR 

NULV 

NWGR 

NWLV 

NZTAGR 

NZTALV 

OLDRHO(NI,NJ) 

OLDZTA(NI,NJ) 

P 

PIM2 

PSI(NI,NJ) 
PSIAUX(NIPA, 
NJM2,JMAXM1 ) 
PSIBFT(NIM2,2) 
PSIM(NPB,NPB) 
PSIO{NJ) 

RDCP 

RHO(NI,NJ) 

RHOO(NJ) 

RI(NI,NJ) 

RIDZSG 

S(NJ) 

SBV 

SPR 

SPSI 

SRHO 

SRTP 
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Maximum number of levels of shading for p graphs. 

Maximum number of levels of shading for log-|g(Ri) graphs. 
Number of time steps. 

Next time step number of which u is to be plotted. 

Maximum number of levels of shading for u graphs. 

Next time step number at which w is to be plotted. 

Maximum number of levels of shading for w graphs. 

Next time step number at which ? is to be plotted. 

Maximum number of levels of shading for c graphs. 

Density from the previous timestep. (old p) 

Vorticity from the previous timestep. (old ?) 

Pressure, (p) 

7 T mul tipi ied by 2. 

Streamfunction . (\fj) 

O 

The solution to V for each vertical level on the 

surface of the varrier, ((ilj^l^)^) 

Fourier transform of the upper boundary of (c^(N^+ 1)) 
Barrier influence coefficient matrix, (a) 

Initial unperturbed streamfunction. (iF°(z)) 

R/Cp = 2/7 
Density, (p) 

Initial unperturbed density. (p“°(z)) 

Richardson number, (Ri) 

CRI / (DELTAZ ** 2 * G) . 

Initial stability. (S°(z)) 

Stability in the constant Brunt-Vaisala frequency case. 
Pressure at the lower boundary, (p^) 

Streamfunction on the lower boundary, (ijj^,) 

Density on the lower boundary, (p^) 

SRHO * STMP / SPR ** RDCP. 


■nil I III ■■■ I in IIIHMIIBI I 
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STMP 

SYM(53) 

T 

TM(NPB,NPB) 

TO 

U(NI,NJ) 

UAl (NJ) through 
UA4(NJ) 
UA5(NJ) 

UA6(NJ) 

UO(NJ) 

U1 

U2 

W(NI,NJ) 

ZETA(NI,NJ) 


AD(NJ) 

B(NJ) 

D(NJ) 

DB 

DT 

DPSIl 

DPSI2 

DRHO(NJ) 

D2PSI(NJ) 


Temperature at the lower boundary. (T^) 

Symbols used in printout of fields. 

Temperature. (T) 

Temporary array formed from PSIM. 

Base temperature for constant lapse rate atmosphere. (T^) 
Horizontal velocity, (u) 

Utility arrays used by subroutines FA, CALCRI, and CALCU. 

Utility array used by subroutines FA and CALCRI. 

Utility array used by subroutine CALCRI. 

Initial horizontal velocity. (u°(z)) 

Horizontal velocity to be added to u at all levels, (u^ ) 
Base horizontal velocity for profiles. (U 2 ) 

Vertical velocity, (w) 

Vorticity. (c) 

Variables Used Uniquely in Subroutine CALCRI 

Diagonal elements of the spline matrix, (a^) 

Elements of the b matrix, (b) 

Backward-difference derivatives of or p with respect 
to z . ( d ) 

Term which incorporates vorticity into bottom boundary 
condition. (a'Az/2) 

Term which incorporates vorticity into top boundary 
condition. (b'Az/2) 

Temporary storage variable for first derivative of aJj. 

Temporary storage variable for first derivative of ip. 

First derivative of p in the vertical. (9p/3z) 

2 2 

Second derivative of ip in the vertical . (3 ip/dz ) 
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R 

UL 

UP(NJ) 

Z(NJ) 


AD(NJ) 

B(NJ) 

D(NJ) 

UL 

Z(NJ) 


FTl (NIM2) 
FTR(NIM2) 
T 


BI(NJ) 

BR(NJ) 

UL 

UP(NJ) 

ZI(NJ) 

ZR(NJ) 


Temporary storage variable for the Richardson number. 
Temporary variable used in inverting the spline matrix. 


Temporary array used in inverting the spline matrix. (Up) 
Temporary array used in inverting the spline matrix, (f) 


Variables Used Uniquely in Subroutine CALCU 

Diagonal elements of the spline matrix, (a^) 

Elements of the b matrix, (b) 

Backward-di fference derivative of ifj with respect to z. (d) 
Temporary variable used in inverting the spline matrix, 
(u,) 

Temporary array used in inverting the spline matrix, (f) 


Variables Used Uniquely in Subroutine CALCW 

Imaginary part of the Fourier transform of t|). (Im(c^(z))) 
Real part of the Fourier transform of ip. (Re(c^(z))) 
Temporary variable used to store Fourier coefficients. 


Variables Used Uniquely in Subroutine FA 

Imaginary part of the elements of the b matrix. (Im(b)) 
Real part of the elements of the b matrix. (Re(b)) 
Temporary variable used in inversion of the relaxation 
matrix, (u^ ) 

Temporary array used in inversion of the relaxation 
matrix. (Up) 

Imaginary part of the temporary array used in inversion 
of the relaxation matrix (Im(f)) 

Real part of the temporary array used in inversion of the 
relaxation matrix. (Re(f)) 
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Variables Used Uniquely in Subroutine FAST 


Except for the addition of the calling parameter INV, this subroutine 
is the fast Fourier transform subroutine from the UCLA BMD library. The 
calling parameters are: 


INV 

M 

N 

X(N) 

Y(N) 


+1 , the Fourier transform is calculated. 

-1 , the inverse Fourier transform is calculated. 

Power of two of N. 

Dimension of the complex array (X,Y). 

Real part of the array to be transformed, and also the 
real part of the transform. 

Imaginary part of the array to be transformed, and also 
the imaginary part of the transform. 


D 

I PI 

N1 

N2 

N2M1 

T 

X(N1 ,N1 ) 

XMAX 

Y 


Variables Used Uniquely in Subroutine GE 

Ratio of the element beneath the pivotal element to the 
pivotal element. 

Lowest row number for which an element beneath the pivotal 
element is to be zeroed out. 

Dimensions of the array X. 

Dimensions of the square matrix for which the determinant 
is calculated. 

N2 minus 1 . 

Temporary variable for the exchange of column elements. 

Array from which the matrix is selected for calculation 
of the determinant. 

Largest element in a row of the matrix. 

Temporary variable for the storage of XMAX. 
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D 

DIFI 

DIFJ 

DL 

DV 

E 

EL 

IDL 

lEL 

I MAX 
IMIN 

LINE2(1 32) 
NEW(132) 

NQ 

0LD(132) 


SIGPRT 
TEMP 
VINC 
VI (13) 

V2(13) 

X{NI,NJ) 

XINC 


Variables Used Uniquely in Subroutine GRAFIC 

Temporary value increment between levels of shading. 

Difference in value from one print character in a line 
to the next. 

Number of levels of shading from one vertical grid point 
to the next. 

LOG! 0(D) . 

Number of print lines from one vertical grid point to the 
next multiplied by XINC. 

Largest absolute value to be plotted. 

LOG! 0(E). 

LOGIO(XINC). 

Base part of the values of the lines separating levels 
of shading. 

Number of levels of shading between zero and XMAX. 

Number of levels of shading between zero and XMIN. 

Logical symbol equivalent of LINE24. 

Minimum-adjusted value at each print character in the 

line corresponding to the present vertical grid point. 

Number corresponding to variable to be plotted. 

Minimum-adjusted value at each print character in the 
line corresponding to the previous vertical grid 
point. 

Significant part of XINC. 

Temporary variable for storage of the symbol index. 

Increment between the values of VI and V2. 

Significant part of the values of the lines separating 
1 evel s of shading . 

Significant part of the values of the lines separating 
1 evel s of shading. 

Array to be plotted. 

Value increment between levels of shading. 
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XMAX 

XMIN 

Y 


BLANK 

BLANK4 

BLANKS 

DASH 

DOT 

LINEl (132) 
LINE2(132) 
LINE4(112) 
LINE5(112) 
LINE6(n 2) 
NAMEl (2,6) 
NAME2(2,6) 
NFB 

SYM1(53) 

SYM2(53) 


None. 


ALPHA 

X(NI.NJ) 


Maximum value in the field; also, maximum value to be 
plotted. 

Minimum value in the field; also, minimum value to be 
plotted. 

Temporary storage variable for X array elements. 

Variables Used Uniquely in Subroutine GRAFIN 

1 blank symbol. 

4 blank symbols. 

8 blank symbols. 

Symbol used in annotation of shading value scale. 

Symbol used in annotation of graph. 

Logical symbol equivalent of LINEl 4. 

Logical symbol equivalent of LINE24. 

Logical symbol equivalent of LINE48. 

Logical symbol equivalent of LINE58. 

Logical symbol equivalent of LINE68. 

Titles of quantities on graphical output. 

Titles of quantities on graphical output. 

First word in LINE24 and LINE34 which is to be filled with 
blank symbols. 

Symbols used in graph of field. 

Symbols used in graph of field. 

Variables Used Uniquely in Subroutine PTSPEC 


Variables Used Uniquely in Subroutine SP 
Barrier influence coefficient, (a) 

Variable for which a superposition solution is required, 
(ij; or p) 
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Value of the array X on the lower boundary, 
or 7^^^) 
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CONSIDERATIONS IN RUNNING THE PROGRAM 


The user must specify all pertinent dimensions at the beginning of the 
first routine. The variables whose dimensions are subject to change are: 


PSI(NI,NJ) 

RH0(NI,NJ) 

ZETA(NI,NJ) 

0LDRH0{NI,NJ) 

OLDZTA(NI,NJ) 

NEWRHO(NI,NJ) 

NEWZTA(NI,NJ) 

U(NI,NJ) 

W(NI,NJ) 

RI{NI,NJ) 

A(NI,NJ) 

PSIAUX(NIPA,NJM2,JMAXM1) 

FTR(NI,NJ) 

FTI(NI,NJ) 

PSIBFT(NIM2, 2) 

UO(NJ) 

S(NJ) 

RHOO(NJ) 

PSIO(NJ) 

where NI, NJ, NPB, NPBP2, NIM2, NJM2, 

in the previous section. 


AD(NIM2) 

AK(NIM2) 

FI (NIM2) 

F2(NIM2) 

AS(NJ) 

UAl(NJ) 

UA2(NJ) 

UA3(NJ) 

UA4(NJ) 

UA5(NJ) 

UA6(NJ) 

BETA(NPB,NPB) 

PSIM(NPB.NPB) 

TM(NPB,NPB) 

LINE34(33,NJSM1 ) 

ITYPE(NI.NJ) 

IB(NPBP2) 

JB(NPBP2) 

LINE3(132,NJSM1 ) 

NIPA, JMAXMl , and NJSMl are defined 
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On the IBM 360/91, the approximate execution space is: 60 * NI * NJ + 

8 * NIPA * NJM2 * JMAXMl + 6 * NI + 11 * NJ + 80K IBM bytes, or roughly, 

8 * (JMAXMl + 8) * NI * NJ + 80K IBM bytes. Approximately 7K bytes of 
additional space are required for each tape used in conjunction with the 
program, for buffering purposes. 

On the IBM 360/91, the approximate execution time is: .00013 * NI * 

NJ * NT * (1 + .0017 * NPB ** 2) +2 seconds. The running time is not 
appreciably increased by moderate increases in the amount of graphical output. 

A version of this program exists which is compatible with CDC systems. 

A SAMPLE CASE 

Consider the Lyra case for u =25 m/sec, T = 250K, and y = 0 on a 

0 0 

64 x 24 grid for 80 timesteps, with Ax = Az = 625 m. At = 10 sec, and a 
2500 m high by 1875 m wide rectangular barrier. Suppose that the entire 
field of ijj is to be output on a vertically exaggerated scale at time step 
60, and ip, p, c, u, w, and log-j^CRi) are to be output at time step 80, with 
no data to be saved. The source program deck, the data card deck, and the 
actual output from this case are displayed on pages 77 through 97. The esti- 
mated execution space and time on the IBM 360/91 are 237K and 17 seconds. 
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Listing of source program 


s _1S THE A I N nOUT I NEJ, 'tiHI CH PLADS THE DATA, 1 N1 T i'AL I Z ES "T HF' 
c FIELOSt and I.NTEGRATEa I HE DEr.SITV AND V J RT tZ 1 T T :£ ■< JA H ^ 

_I IMJ^LICIT REAL*8( A-H.O-2» 

“REA1_>8 P8I( 66.24).WHJ( 66 • 24 ) t 2ETA( t6.24lV " 

1 CLDRHO( 66 ,24> ,yLDZTA( 56 , 24 ) , NEWRmO { 6bt24i. 

2 NEWZTA( 66,24),F2t 64). 

.2 U< 66,2*J,W( 55 .24 > . R I ( _6e .24} , A( 66.24). 

3 PSIAUXC 67 ,22. * ) . = TS < 65.24).FTI( 6o . 2 4 J , ^ a 1 3F r't ' 6 4. 2 > . ~ 

4 UO(24) .S( 24) , ?HJO( 24) .PSI 0 t 241 . AO( 6 a).A<( 64).F1( 64), 

5 AS<24) ,UA1 (24).JA3(24),JA3(24)*UA4(24).JA5(24).UA6(24). 

_ 6 BFTA{ 10.10 >.PSlM{t0.10) .TH(IO, 10), 

7 L INE43 ( I A ) ,HnC 58 ( 1 4 ) . i. INR68{ 1 4 I 

REAL* 4 LINEI4(33),LINC2 4<33 ) . L t NE34 ( 33 , 2 ) , 'aA4E(2.b} , 

I PLANK4 

*J:<J?GEfi*» IPS IGR t 20 ), IPH3GR(20 ), IZTAGP { 20 ), IJ6R( 20) . I »GR{ 20) , 

T ■■ ■ “I PI (20) ;7TVPEl “6B',’4T . 1 B( 1 2) , J3( I 2) , iOJ Mi e“) 

LOGICAL*! L1NE3C I 32.2) ,SYM( 53) 

EQUIVALENCE ( L I NE 3 ( I ) . L I NE 3 4 I 1 ) ) , (NEiKZTAI l).J(l).«( 1).RI( 1)), 

I __ _ _ ( PS! (1 ). = IRC 1 )),( NEWRHQI 1 ) .FT I ( 1 ), Al 1 ) ) 

C OMMON/Cl/OEL TAX,DXM2.0ELTAZ.02M2.DZ5Q,OZD6.G,Ni4l',NJ4l i 1 EXP2 , 

1 I S. JS, I E. JE‘. ISit • lEW .flPDPl ,JMAX, ISO, iEO 

C OMMQN/C 2/LI NE 4 8, LI NE 58, LINE 68. LINE 14 ,LINC2 4. NAME , 

1 MT. IN.JN. 1N81.N1S.M SM1 .1L,SYM 

C SPECIFY PHVSI CAL factors ■ 

G = 9,8 

CMTMPS = I8S2, / 3600. 

SPSI = 0. 

SRHO = ■|.*25~~ 

ST MR = 2 73. 

S3R = 1 . OD+05 

SRTP = SRHO ♦ STMP / *♦ <2./7.) 

OAlR ■= ■9,'76T5^03 — — 

RDCP = 2. / 7. 

c 

C A_, RE^, _y»RIIE.SPACE_AND _T IME_(>ATA 

RE AD (5. 1001 ) NI.NJ.NT.N^a.IBT.I SAVE.OELTAX.DELTAZ.OE^TAT 
10)1 FORM AT { 6 14 . 1 P3DI ) ,2 1 

*> R I TE ( 6 . 200 1 ) N1 , N J ,M , NPO , I BT , 1 S AVE, DELT AX , OELT AZ. JELT AT 
20)1 “FORMAT ( 1 HI ,• summary )f iM^jt DATA FOR NON -1 T E R A T 1 V E . SOLID «■. 

1 'PARPIER oRUGRAM* / /// • Nl =',I4,«, RJ =*,i4,", NT := • , I 4 , 

2 *, N^3 =*.I4.", IDT =«,I4.». ISAVd =*,1A,«, OX =», 

3 10010.2, •, DZ = • , 1 PDl 0. 2. ■ , DT =*,1»010.2) 

C CALCULATE QUANTITIES H I C H DEPEND ON SPACE A>1U IlME DATA 

NIMl = NI - 1 

N IM2 = NT - 2 

Z CALCULATE POWER OF TWO OF NUMBER OF UNIQUE HJRIZJNlA- GRID POINTS 

I EXP2 =■ 1 " " “ 

1 = NIM2 

5 1=1X2 

I F ( I . EQ . t ) GO TO 10 
IEXP2 = ieXP2 41 ^ 

SO TO 5 

10 NJMI = NJ - I 

NJM2 = NJ - 2 

MI = IBT 
lET = IBT ♦ NT 

1B2S = O 

NP3M1 = NPB - 1 
N^UPl = NPB + 1 

NPUP2 = NPB 4 2 

DXM2 .=. DELTAX._*._2 . 

DXSQ = OELTAX 4* 2 
0ZM2 = DELTAZ * 2, 

DZD2 = DELTAZ X 2. 

DZSQ = DELTAZ,. ♦.*..2 _ 

DZ06 = DELTAZ / 6. 

DZ2DX2 = 2. * OZSQ X DXSQ 

FAC = -I, X 112. 4 DELTAX 4 OLLTAZ) 

. _ - FACM2„=,.FAC . 4 .. 2., _ 

FACM4 = FAC 4 4, 

FACM43 = FACM4 X 3. 

DELTAT = DELTAT 4 2. 

.c 

C 8. READ, IkRITE BARRIER DATA 

REAO( 5, 1002) ( H)( I I . JU( 1 ) ,1=2.NP3P1 ) 

10)2 FOR,MArt20I4J . . . 


MR I TE (6.2002 } ( M((l. JB(I) .1-2.NPBP1) 

2002 FO^MAK IHO. 'COORO INATES OF tARRIER . 1 0 ( 2X . • ( * • I 2 . • . • . 12 , • I • | | 

C CALCULATE OUANIITIES AHICH DEPEND DN BARRIER DATA 

I WD = IBtNPBPlI - lB t2l 

NIPA = NIM2 ♦ IWB " 

C CALCULATE TYPE OF EACH GRID POINT 

CALL PTSPECd TTPE.IB. J0.NI,NJ,NPBP21 

JMAXHI = UMAX - t 

NOUH = 2 * HODCNJ A( (N 1 *7 ) /2 4-1 I '•’N [ M2*3 4-NlP AAM JM2» Ji4 A Xm -t-NPB**? • A ) 
DO 15 I = l.NDUM 
15 IOUM( I> * 0 

Z C. REAdT” W^Te^GR APHI CAL OUTPUT DATA 

C 

READ! 5, 10031 IS.US.IN.JN.NJSMI . NPS I L V .NRHOu W' , N^ T Al V . NUL V , NWL V . 

-I NRlLJ/.JP_bJGRj_IiiHUGR. J2JAGP.IUGR, 1*l.R. 1 HI GR 

1003 FORM AT( 514/61 4,6< /20I 4 1 ) 

liRITE<6.20 03> lS.JS,IN,JM.NJSMl,NPSILV.NRHaLV.NZrA_V.MULV .NWLV, 

1 NRI LV.I PSIUR . IRHOGR. IZT AGR, lUGR , 1* GR. IR IGR 

^20JJ -FXIHMAT_|JHO. J4. <5 . = J 4 , -JJM _J,V 

1 •• NJSMl =*,14 //• NPSILV =».14.*, NRH3i_v =".I4, 

2 •• N2TALV »‘.I4,*, NULV =*,14,*, N*LV =*.14. 

3 *. NR ILV =»,I4//» IPSIGR =',20I4//» IkHJGH =*,2014// 

A J_a2TAiB_s_*.20J.4/21! — LUGR. !^0-l4yy^ »_4 * jR..= t+JSOlA-/-/- 

S • IRIGR = '.ZOlAl 

C CALCULATE QUANTITIES WHICH DEPEND ON GRAPHICAL DATA 

NPSIGR = 1 

NRHOGR = 1 

NZTAGR =1 
NUGR = 1 
NWGR = 1 

NRIGR_=_1 

1£ = IS ♦ IN - 1 
JE = JS ♦ JN - 1 

C calculate HCiHIZCNTAL EXTENT OF GRID POINTS AT lirtICH W IS TO BE 

C CALCULATED 

ISW = IS 

IF< IS.EO.I > ISW = 2 

1 ew JE 

IFIIE.EO.NI) lEN = NIMt 
I NMl 5= I N - 1 

Z SET UP ANNOTATION FOR GRAPHICAL OUTPUT 

CALL GRAFINIL I NC3 4 . L I NE3 . N JS M 1 ) 

IF(iaT.NE.O) GO TO 230 

c 

C O. READ WIND VELOCITY DATA. AND CALCULATE UO 


READ(5.1004) iCASE 

1004 FORMAT! 14) 

GO TO ( 60 . 70, 30. go, I 0 0, 1 1 0 » . I CASE 

C SOUNDING CASE (CALCU-AT E RHf .O IN ADD I TIO N T3 UOl 

60 fcRI TE C6 .2004 ) 

2004 FORMA T( 1 HO. ■ SOUND I NG DATA C ASE » //6X , • U • . I 1 X . • T • , I I X , • P • / > 

DO 63 J = 1 ,NJ 

READ(S.i005) UOim.T.P 

1005 F0RM4T< 1P3D10.2) 

WRI1E(6.2005) U0(J).T.P 

20 0 5 F ORMATI 1X,1PD10.2.1P2 012.2) 

___J0(JJ - U01JJ_.*_CIUMPS - . - - - 

T = T ♦ 273. 

P = P ♦ I .ODf 02 

65 RHOO(J) = SRTP * P ** (2./7.) / T 

GO TO 12 0 

z constant velocity case 
70 READ! 5. I 006) Ul 

1006 F OHM A T ( 1 PDI 0. 2 ) 

WR I TE( 6, 2006) _Ul 

2006 form ATI IHO. 'CONST ANT VELOCITY CASE, Ul =».IPDU).2) 

DO 75 J = 1 ,NJ 

75 U0( J ) = Ul 

GO TO 120 

C CONSTANT SHEAR CASE 

80 READ( S, 1007) Ul ,C 

1007 FURMATf IP2D10 .2 ) 

ir.'RI TE (6,2007) Ul.C . _ . 

2007 FORM\T( IHO, 'CONSTANT SHEAR CASE, Ul = • ,1 PJ 1 0 . 2 . • , C =*,IPD10.2// 

1 6X. 'U' /) 

COZ = C * DELTAZ 

DO 85 U = l.NJ 

U0(J) = Ul 4^ CTZ ♦ IJ - 1) 

85 aR ITE(6,2008) U0(J) 

2008 F JRM AT ( I X. IPD 10 .2 ) 

Ga. IO .12.0 — - . 
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C EXPONENTIAL PROFILE CASE 

90 REAOtS.lOOS) U1.U2.C 
WRITE <6 .2009) UI.U2.C 

2009 FJR M AT! I MO , • EXPO NENT! A L PR QFJj.^ C A SE « U1 _=• * J P31 Q , 2 . « . U ? = « . _ 

I lPDtO.2.*. C =•. IP010.2//'6X, "U*/! 

COZ = C * DELTA2 
00 95 J = T.MJ 

O0{ 3 L_=_JL»1 » U2 ♦ DEXPtCO Z» lJ-l I I 

95 «RI TE(6. 2008) UO(J) 

SO TO 120 

C HYPERBOLIC TANGENTIAL PROFILE CASE 

100 READ(5.100a) U1.U2.LI.M3Z 
lOOa FORMAT! IP2O10. 2.214) 

hrRI TE (6 .201 O) U1.U2.LI.M0Z 

_^OX-a_E-aRMATl-l-HO.._!M YP.Ea ROL 1 L TANGENIJAJ PROflLE-XA SE. jJJ. _=i .JPJ31 Q.2,. 

1 •. U2 =• . IPOl 0.2. * . LI =*.I4.». MOZ =• . 1 4/ZoX. • U* /) 

DO 105 J = I . NJ 

UO(J) = U1 ♦ U2 ♦ DTANH! OFlQATI IJ-LI ) ^MDZ ) ) 

,_ia5 WRITE 16. 20.05JL_UflUi 

GO TO 120 

C OTHER VELOCITY CASE 

I 1 0 WRITE (6.201 1 ) 

__2.01J FORMAT! IMP. «S PEC1F1EJ VELOCITY C A SE » //6X. » J* / i 

REAOI5.1009) <U0( J).Jsl.MJ) 

1009 = ORM AT ( I P8O10 .2) 

OO 115 J = l.NJ 

l_13._WHJI£i6. 20 jQBJ UO J_J_) 

c E. calculate PSIO 

c 

. -J20 P5I0! 1 ) = 0. 

DO 125 J = 2.NJ 

125 PSIOIJ) = PSIOIJ-1) - 02D2 ♦ (UOIJ-1) ♦ UO!J)l 
I F! ICASE.EQ.t ) GO TO 170 


C F. READ. WRITE STABIi-ITV DATA, AND CALCULATE STAJI_ITY 

READ(S.1004) JCASE 

GO TO 1 i.3D. 1A0 j_15Q_LlJXA5E 

C CONSTANT LAPSE RATE CASE 

130 READIS.1007) GAMMA, TO 
fc R I T E !6 ,201 2 ) GAMMA, TO 

_2012 FORMAT! IHO, »crN,STA NT L APSE RA TE CA5 F, GA M.M A_ ^ , 1 P J 1 0 ^2 •_* 3 .*, 

I IPOIO.2) 

OMGOTO = IDALR - CAMMA) / TO 
DO 135 J = 1 , NJ 

1.3 5 SIJ) ^ DMGDT^ 

GO TCJ f60“ 

C CONSTANT 9RUNT-VA I SALA FREQUENCY CASE 

140 REAO(5,1006) BV 

_ _>R1TE!6.2013J_ QV 

2013 FORMAT! IHO, ‘CONSTANT BRUNT-V A I SAL A FREQUENCY CASE, BV = • . I PD 1 3 . 2 ) 
SBV = BV ♦* 2 / G 

DO 145 J = l.NJ 

i4S„si JJ ^_sa« 

GO TO 160 

C constant RICHAROSCN number CASE 

150 REAOI5.1006) CR 1 

hRlTE!6,20l.4 )_CBI 

2014 FORMAT! IHO. ‘CONSTANT RICMARDSON NUMBER CASE. R1 =*.IPD10.2) 

RIOZSG = CRI / (OZSQ • G) 

S(l) = RIOZSG 4 !U0!2) - UO ! 1 ) ) *4 2 

DO _L5 5_.J_^2i£IJf.’ll 

155 SlJ) = RIOZSG * (UOlJ+1) - UOIJ-1)) 4* 2 / 4. 

SINJ) = RIOZSG 4 lUOlNJ) - UOINJMl)) 44 2 
C 

C G,._CALCULjAt£_BHaQ 

160 RHOO! 1 ) = SRHO 

DO 165 J = 2.NJ 

165 RHQOlJ) = BHOOIJ-11 4 D£ I - DZ 02 4 ! S ( J- 1 ) ♦ S ! J ) i ) 

C INITIALIZE VQRTICITY. ANGU.AR ARGUMENTS, MATRIX DIAGONAL 

C elements. ANO FOURIER TRANSFORM PSl = O CN TOP 3JJNJARY FOR 

C_ CALCULATION_£lF_£ARRI£iT_I>.FLUENCE- coefficients ... . 

170 DO 17S I = l.Nl 

JO 175 J = l.NJ 
175 ZETAt I, J) = 0. 

P1W2 . = ,_6 .2831 853 07-180 4 00 

DO 130 I = 1.NIM2 

ARGl = PIM2 4 II - 11 / NIM2 

ADI II = 2. 4 DZ20X2 4 II. - OCQSIARGI)) 

. AKII) .= DSINI ARGI-J-/. OELTAX ._ - 


i I I .It UlMJ I U I , II II : UU>J 


00 180 J = 1*2 
1 BO P SIBFT( I . Jl SO. 

C OtFIME SPLINz JilATRlX DIASONAL ELEMENTS 

ASI-I 

DO i82 J s 2.NJM1 
182 AS(J) = A. 

AS(NJ) s 2. 

TT’.^s'olvc POISSON EQUATION fciTii unT t voRfi'cTT 7 'JiTt'j 5b aVc e ea^h 

VERTICAL LEVEL FOR WHICH THERE IS A RDiNT UN THE SJH-ACE OF THE 
BARRIER. AND STORE RESULTS INTO PSIAUX 


DU 190 K = 1. JMAXMl 
ZETAl IB< NPBPl I .K+ 1 I = 1. 

CALL FAIPSI.ZETA.FTR.-TI.PSIDFT.AO.UAI .UA2.UA3 .UA«.UAS.NI.NJ,NIM2, 

... 1 NJM2J 

ZETA( 1 B { NPBPl ) . K+ 1 1 = 0. 

DD 187 J = 1. NJM2 
DO 185 1 = 1 . NIM2 

.185 PSl A JXil..JjXJ_s. P-SJ.ll.JtlJ 

DU 187 L = 1 • IWB 

187 PS I AUXT NIM2+L . J.K I = PS I AUX <L . J . K I 
DO 189 J = 1 . NPB 

_IF< JB( Jj-1 1.NF.K4-1 1 GD TD 1B9 

L = IB(NPBPl) - IEJ(J+1) 

DO 188 I = 1 . NP8 

PSIMd.Jl s PSK IB( IFI )«_. JU(14.1 ) ) 

ISB.TMd.JJ =_PS1MU_.JD .. 

189 continue 

190 CONTINUE 

calculate determinant of alpha matrix 
..CALL GEiIM.DEJL.N.RJl.N.R±L.J.AXCiJ 

1. CALCULATE BETA 

DO .205 X..= .1.. NPB ... _. 

DO 205 J = 1 . NPB 
CFAC = C-1) *♦ ilFJ) / DET 
SELECT COFACTOR MATRIX ELEMENTS 

do” 2’o'o~K~^”i7’nP8 

IFII. EQ.K) GO TO 200 

M = M ♦ 1 

DO 195 L - 1 . NPB 

IFIJ. EO.LI GO TO 195 

N = N ♦ 1 

T MI M.N) = PSl MIK.Ll 
195 continue 
200 CONTINUE 

_ __C ALCULATE _DF.r£RMI.NANr_DF COFACTOR__MATRiX 

C ALL GEITM.CFJET.NPG.NJB'll .IACC2) 

205 3ETA(J,I) = CFAC * CFUET ♦ 10. I I I A CC2- 1 A C C 1 1 • 50 1 

CALCULATE FOURIER TRANS = 09M OF PS I ON TOP BOUNDARY 

DO 20 9 I = ^1,_N1M2 

209 PSIBFTI I .1 ) = PSIOINJl 

CALL FA STI PSIBFTI I .11 .PSIBFTI 1 .2) . N IM2 . I E XP2 . 1 i 

J. INITIALIZ E PSt . R H Dj_A NILXEXA 

DO 210 I = 1 . NI 
DO 210 J = l.NJ 

...21 0 .^SII I .J) = PS lOI J > 

C SUPERPOSITION solution FOR PSl FOP TIME STE=» 0 

CALL SR I P5I .P5I AJ X , IT YPE . BET A . IB.JB,5PSI.NI.NJ.N1M2.MJM2.NIPA, 

1 JMAXMl , NPB, NRBP2) 

_ 0 0 215 1 = l_iJJ I 

DO 215 J = 1 , N J 
DD 214 K = 2.NI 

214 1 FI D ABS IPSI I 1 , JM .LE.DADS I RSIOI K I » ) GO TO 215 

215RH0(I.J) E_(f>SI(I,JJ PSIO IN-I.J J * I RhCO ( <i_. - R i jJ I <- 1 J ). y 

I {PSIOIKl - PSlO(K-l)) + PHQO(X-l) 

C CALCULATE ZET A FROM PSIO 

OJ 220 J = 2,NJM1 

220 ZETAIl.J) = IP5IJIJ -11 t ? S10IJ»M - 2. * P S I 0 1 J JJ / .U.Z S.D 

ZET A( I . 1 ) = ZETAI 1.2) 

ZETAIl.NJ) = ZETAI I .NJMll 
DO 225 J = l.NJ 

DO 225 I =_2 . NI 

225 ZETAIl.J) = ZETAIl.J) 

C ZERO CUT VOPTICITY IsSIOE EAf.RIER 

DO 227 I = ISB.1E3 

DO 22 7_ J = XM1___^ 
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227 I F( irvPEf l»J) .E0.9I ZETA(I,J) > O. 

GO TO 23S 

K • RE AO . PS I.,_RHO.,_^ELAjLjil±0-Qj_PSJLA.U_)^_B£TAa ..lJJf:^j_.A3 .-AND. PSI 8FT 

FROM TAPE OR DISK 

230 REAO(l) PSI .RHO. ZETA. RHUO .PSI AUX.ElETA • ITYPE . AO .PSI JFT ■ 

1 < tom ( I > . 1 = 1 . s ajMi 

2 35 IF( IB25".EQ.I ) GO TO 2 80 

L, CALCULATE GRAPHICAL OUTPUT FOR TIME STEP MT 


1 F(MT .NE. IPS I GR( NP5IGR > I GO TO 24 0 

CALL GRAFICIPSI .A.LINEOA.LINEB.l.NPSILVtNI.HJ.NJaHli 
N^SIGR = NPSIGR ♦ 1 

2A0_ IFIMT ,NF. IRhOGRINRHDGR) ) GO TO 245 

CALL GRAFIC(RH0.A.LlNE34.LlNE3«2.NRH0LVtM.NJiNj9Mli 
NHHOGR = NRHOGR I 

245 I F( MT .NF . IZT AGR (HZTAGK » I GO TO 250 

-CALL GRAF1CIZET A.A, L I\E3».LINE3.3 . NZT A L V. N I . N J . N J a >4 1 I 
NZTAGR = NZTAGR *■ I 

250 I F( MT ,NE. lUGR (NUGR) ) GO TO 255 

CALL CALCU(PSI.ZETA.UiAS.UA1.UA2.UA3.0A4.NI.NJ| 
call GRAFICIU.A.L INE34 .L I NE3.4 .NJLV.N i.NJ.NJSMlI 
NUGR = NUGR * I 

255 I FI MT .NE. IHGP INWGR J ) GC TO 260 

_ CALL CAI CW(PSl.\lt.Fl .F?.AK1. M.NJ.N 1M?1 

C all GR AF IC( W .a ,L 1NL34.LI NE3i5iNhLV.NI .NJ .NJ5M1) 

N1.GR = NWGR ♦ 1 

260 I F t MT .NF. I RI GR (NR I GR) ) GO TO 265 

- .call CALCRKPSI ,RH0fZtTA,RI,AS.UAl,UA?jUA3,JA4,JA5,JA9.NI,NJ) 
CALL GRAFIC{RI.A,LINt34,LlNF3.6.NRILV.NI.NJ,.JSMI> 

NRIGR = NRIGR + 1 

265 I F( MT .L T . lETI GO TO 270 


M, MHITE PSI, RH3, ZETA. RHOO, PSIAUX, BETA, ITTPE. AD. AND PSIBFT 
ONTO TAPE OR DI S< 

IFT I SAVE.EQ.IJ WRI TEI 2 1 ->^3 I j RJlO. Z£T A , R FOO , P5 J A JXj BE T A , 1 T YPE. AO. 

I RSIBFT .( IDUMI 1 ) ,1=1 ,NOjM) 

STOP 

270 I F I MODI MT.25) .NE. 0) GO TO 280 

DEL tat. --DCL.LA T__Z_4. 

N« SET OLDRHO = RHO, ULDZTA = ZETA 

DO .2 7 5 _L_= . 

DO 27S J = I.NJ 
DLDRHOII.J) = RHOd.JI 
275 CLOZTA(I,Jl = ZETAII.Jl 


O, CALCULATE JACOBIANS JIZETA.PSI), JIRHO.PSl) 

230 DO 345 I = 2.N1M1 

D3_ 34 0 3 - I .NJ 

k = 1 TYPE! I. J I 

C NEWRMO FUR POINT TYPES I THROUGH 9 

NEWRHOIl.JI = OL3RHO(I,JJ 

IFIK.NE.O )__GO _T_p _2 85 

C JACOB IANS FOR POINT TYPE 0 

A I = PS I I I-I . JF I ) 

A2 = PS I I I, JF I 1 

AJ_ = _PSI ( IFJ , JAJ_I 

A4 = PS I ( I-l . J1 
A6 = PS I < I FI . J| 

A7 = PSII I-l, J-l» 

A 8 . = . PS I .( I , J-.3.I 

A9 = PS I ( I Fl , J-1 ) 

Bl = A8 F A9 - A2 - A3 

B2 = - A7 - A3 F A1 F A2 

B3 =_A6 F_ A3 - A4 - AI 

B4 = - A9 - A6 F A7 F A4 

B5 = A6 - A2 

36 = - A8 F A4 

P7_j=_A 2 ri_A4 

S3 = - A6 F A3 

AJZP = FAC F (31 ♦ Z£TA(IF1,J) F B2 ♦ ZETAIl-l.Jl F 33 ♦ 

I ZETAd.J+Il F 34 F ZtTA(I,J-l) F ri5 F E^TAlIFl.JFll F B6 4 

.2 2ETA( I - I > 3^1 J—F .E7._*,.2ETA( I - I . J F 1 ) F j 8 * ZLTA ( I F 1 , J-IJJ 

DRHO = (RHO(IFl,J) - RHO(I-I,J)> X 2. 

C CALCULATE. NEwRMO FOR POINT TYPE 0 

NEWRHOII.JJ = FAC F (Bl F RHOIIH.J) F 82 F RHJ(l-I,J) F B3 * 
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1 RH3fI.j4-l) ♦ B4 » RHO(I.J-ll «■ B5 * AH3 < 1 4^ 1 » J*- 1 } * 

2 B6 ♦ RHClI-I.J-ll + B7 * RHO({-it J41 J ♦ 88 • 

3 RH3 (141 , J-1 ) ) * DELTA! 4 OLOHHDd.Jl 

._GQ..JI1J33S 

JACOBIAN FOR POINT TTPE 1 
28S IF(K,NE.1J GO TO 290 
A5 = PSKI.J) 

AZ.._= PSI-1 1,-1 

AS - PSKI.J-rl) 

A9 = PS I ( ( 41 , J-l 1 

AJZP = FACH2 ♦ (ZFTA(I41,J) ♦ ( A8 4 A9 - 2. ♦ A5) 4 ZETA(I-I.JI * 

1 lri._A7..^-A3 _l-.a-..-«-ASJ. 4.-ZE.TA1I . J-l J._* _X-_A9.4 _A7J > 

2 < ZET A< I-l . J-l 1 - ZETA( 141 . J-l) 1 * (- a 8 4 Aal) 

DRHO = Om 
GO TO 335 

JACGDIAN.FQR.£DUtU_rXE»E_2 

290 1F(K.NE.2> CO TO 295 
A1 = PST(1-1»J4I) 

A 2 = PSI (I .J4l) 

A3_ = _PS J 1 1 41 ^-4_U 

AJZP = 4ACM2 ♦ (ZETA(l4l.J> * < - A2 - A3) 4 ZETAd-l.J) ♦ 

1 (Al 4 A2) 4 ZETA(lpJ41) ♦ CA3 - All 4 A2 4 ( ZE 1 A ( 1 - 1 « J 4 1 ) - 

2 ZETA(|4I,J41)|) 

DRH Q = O. 

GO TO 335 
295 K M2 s K — 2 

GO TO < 300,305.31 0,31 5. 320*325t330) .KM2 

_ .JACOeiAN FOR POI NT T T PE 3 

300 A 1 = PSI ( 1-1 . J41 ) 

A4 = PS I ( 1-1 , J) 

A 7 = PSI ( I-l , J-l ) 

A JZP. = F ACM2L *_<Z£TALJtJ ..JJ__A__1-_AZ_4_ AI ) _t ZL T4 1 J j J4LJ j» _ 

1 <- AA - Al ) 4 ZETA(I.J-l) * <A7 4 A4 ) 4 A4 • 

2 < ZETAl I-l , J-l ) - ZETA( 1-1 . J4 1 } 1 1 

DSHD = PHOd.J) - «HQ(I-l,J) 

..GO T0_335 

JACOBIAN FOR POINT TY^E 4 
305 A3 * PS I ( I4I , J4I I 
A6 = PSI (141, J) 

A 9 =..P3I II.4J, Jrll , . 

AJZP = FACM2 » (ZETA(I41.J) * ( A9 - A3» 4 ZErA(i.J4l> * (A6 4 A3) 

1 4 ZETA(I,J-l) * (- A9 - A6) 4 A6 ♦ ( ZE 1 A ( 1 41 . J 4 1 } - 

2 ZETA( I 41 , J-l )) ) 

__ .DRHCL- ,-RriQ(l4j..J ) -RH01J-.-U 

GO TO 335 

JACOBIAN FOR POINT TYPE 5 

310 AJZP = FACMA * PSId-I.J + l) * (ZETAd-l,J) - ZErA(l.J41)) 

3arto_s:_0, — 

GO TO 335 

JACOBIAN FOR POINT TYPE 6 

315 AJZP = FACM4 * PSI(141.J41) * (Z£TAd»J41) - Z£TAd4l.JI) 

DRH0_=_0« 

GO TO 335 

JACOBIAN FDR POINT TYPE 7 
32 0 A 1 = PS I d-1 , J4 1 ) 

A2_= .PSd I ,J4JJ 

A3 = PSI (141 , J41 ) 

A4 = PS I ( l-l , J) 

A 7 = PSI ( I-l , J-l ) 

AJZP = FACM43 ♦ (ZETA(141,J) ♦ (- A2 - A3) 4 ,.oTA(l-l,J) • 

1 (- A7 4 Al 4 A2i + ZETA(I,J41) * (A3 - A4 - Al ) 4 

2 ZETAd.J-1) 4 IA7 4 A4 ) 4 A2 4 < ZE TA d - 1 , J4 I ) - 

3 ZET A ( 1 +1 ,. J 4JJJ_4_ A4 4 (ZETAd-l . Jjr^l J_t_ZE T A ( I - 1 , J 4 1 )J J 

DRHQ = (RHO(I41,J) - RrlO<I-I,J)l /• 2. 

GO TO 335 

JACOBIAN FOR >^0INT TYPE 8 

325 A L^PS I ( , J4 1 1 

A 2 = PS I ( I , J4 i ) 

A3 = PSI ( I4l , J41 ) 

A6 = P5Id41,J) 

A 9 =.^PS I (.I41..J-1L. . ..... 

AJZP = FACM43 4 (ZETA(IM,J) 4 ( A9 - A2 - AJ) 4 ZErA(l-l.J) * 

X (Al 4 A2» 4 ZETa{I,J 41) 4 (A6 4 A3 - Al) 4 ZElAd.J-l) 4 

2 <- A9 - A6 ) 4 AS 4 ( ZET A ( I 4 1 , J 4 1 > - 2£ T A ( 1 4 1 , J - 1 ) ) 4 A2 4 

^3 ^lZETA(,(-_lj J4.lJ__.r_ ZETA( I 41 . J+IJ ).J . „ 

DRHO = (RMO{I41,J) - HH0(I-1,J))/ 2. 

GO TO 335 

NEWZTA FOR POINT TYPE 9 

330 NEV.'ZTA{ 1 . J) = O LDZTA(l. J) 

GO TO 340 

CALCULATE NEWZTA FUR POINT TYPES 0 THROUGH 8 
335 NEWZTA(I,J) = (AJZP - (G 4 OHHD) / (DELTAX 4 N.IJJ(J))) 4 DELTAT 

I .+_oi_:>zT.A(.i_.ja 
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340 CONTINUf; 

3*5 CONTINUE 

1 F(MOO( MT.25) .EQ. O) SO TO 355 
P“.~SeT OLOZTA = ZETA* OLD^HO = RHO 
DO 350 I s Z. NIMl 

DO_350 J =_l_tJNJ 

a'LDZTAil.jl = ZCTA(I.J) 

350 OLDRHO(I.J) = RHOd.J) 

GO TO 365 

355 DELTAT = PELT AT » 2. 

I F< I825.E0.1 » GO TO 300 
I 825 = 1 
MT = MT - I 

. __iJO_T9_J65 

360 1B2S - O 

O. SET ZETA = NEWZTA. RHO = NEtfRHO 

3^ 5 30 Xr O 1 s 2, NI Ml ' 

DO 370 J = l.NJ 
ZETA(I.J) = NEWZTA(I.J) 

j.70_PHO(i J jj 

DO 375 J = 1 . NJ 
ZCTA(1«JI = ZETAtMMl.JI 
ZETA(Nt.J) = ZETA(2,J) 

WHO I ! I _=_BiH_0( Nl Ml ,J ) 

375 RH0(NT.J|‘= R)jOl2.J) 

R. CALCULATE PS I FOR NEXT TIME STEP 

SOL’vIi^P trrs^N~Fo U ATIUN FDR~'PSI ’ 

CALL FA(PSI.ZETA.FTR,Fri,P51BFT,A0,UAl,UA2,JA3.UA4,UA5.NI,NJ,NIM?, 
1 NJM2) 

SUPERPOSITION solution FOR PSl AT TIME STEP MT 

C ALL SPtPSI.PSIAUX.ITyPE.BETA.IB.JB.SPSI.NI.^J.MlMZ.NJMP.NIPA, 

I JMAXMI ,nP6.NpBP2) 

MT = _MT_*_1 

GO TO 235 



|SJnPQUTlN£ CALCRII IPSI »RH0tZETA,RI.An.UP«RfZ,D.D2PSl tDRHD.NI ,NJ) 

' This RP'uflNl TALrULrfrX~I Hr~Rl CHXPDS'Cr^ KUMi3ER~~FTZi-J Hi ' 

IMPLICIT REAL4H1 A-H,D-Z» 

_ __?£AL*3 P5I< NI ,NJ) ,WH0( NJ ,NJ) jZETA(NI, NJ), RI {Nl ,MJ> • JPl NJI j 

t AD(NJ).8(NJ),Z<.NJ).D(NJ),D2PSI(NJ).0RJ0{NJ) 

CQMMON/Cl/OELTAX, DXM2 ♦ D E LT AZ . D ZM 2 , O ZS U . DZ 06 . 3 . N I * I , ■» J M I . I E X»2 , 

I IS.JS.tE.JE«IS«.lEw>NPBPl.JMAX,lsa.lEO 


A. CALCULATE SECOND DERIVATIVE OF PSl 
UP( 1 ) = ADIl » 

p_( LL = _ Q . 

DO 10 I = IS. IE 

boundary CUNDITIQNS on SECOND DERIVATIVE OF ^Si 
D2P5I 111 = ZETA( 1.1) 

D2PSKNJ) = ZCJ_AJ.I_.NJJ 

OS = D2PS 1 ( 1 ) * D7D6 

DT = D2PSHNJ) * OZ06 
CALCULATE 8 MATRIX E-EMENTS 

D0_1..J..= ..2.NJ 

D{J) = <PSI(I«J) - PSl(I.J-ll) / oeltaz 

1 BlJ-l) = 3. * (0(J» ♦ D(J-l)) 

QtlJ = BID - 3. ♦ OB 

8(NJ) = 3. ♦ tD(NJ) * P TI 

TRIDIAGONAL SPLINE MATRIX INVERSION SCHEME FjM = 1 RST DERIVATIVE OF 
PSl 

Z ( D = B( 1 » 

DO 2 J^2,NJ 

JL = -1 . / UP ( J-1 > 

U^(J) = AO(J) + UL 

2 Z<J) = B(J) + UL ♦ ZIJ-l) 

DPSIl =. ZINJ) ^._UP(NJ) _ 

CALCULATE SECOND DERIVATIVE OF PSl FROM FIRST DERil/ATIVE OF PS I 
DO 3 J 3: 2.NJMI 
K = NJ-J+l 

DPSI2 = (Z{K ) , r_DP.SI 1 »_>f_UP<.KJ 

D2PSKK) = (-4. ♦ OPSI2 F 6. ♦ D{KF1) - 2. • PPSID / OELTAZ 

3 DPSI 1 = 0PSI2 

B... C A LC ULAI.e__EJaST_ iDEi LVA.T_LyE_X)l^_RHD 
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Oa 4 J s 2.NJ 

3(J» * IRHO(I.J| - R-iO(I.J-l)| / OELTAZ 

Z C ALCi;LAl£J_^TiUJL_£:LiMf.?irS 

A 3(J-l> = 3. ♦ (r>(J) f OlJ-lti 
B(NJ) = 3. * D(NJ) 

; TRIOl AGONAL SPLINE MATRIX INVERSION SCHEME 

I II I = BII I 

OO 5 J = 2.NJ 
UL = -1. / U^( J-1 ) 

JP(J) = AO(J> ♦ UL 

5 21J> = BIJ) ♦ UL * ZO-I) 

DRMOINJj = ZINJI / U^INJl 
00 6 J=2.NJ 

K = NJ ._- _1 

6 DRHOIK) = (2IK) - DRHDIK^-I)) / UP ( K ) 

C. CALCULATE RI FROM ORH 0 . D2PSI 


DO 9 J = JStJE 
1FID2PSI r J) .EQ.O. 1 GO TO B 

R = - G * DRHO(J) / tRHD(ItJ) * 02PSI(J)**2) 

IFIR. GI^ -25 ) GO TO 7- 

R = . 25 
GO TO 9 

7 IFIR. LT. 9. 99) OO TO 9 

8_a,_=_S.9:9 

9 R I ( I . J ) = DL3GI0I R) 

1 0 C ONT I NUE 
RETURN 

ENO ^ 

1SU3R3UTINE CALCUK PS I . ZET A , J . A 0 . OP . B . 2 . D »N I , NJ ) 

THIS ROUTINE CALCULATES THE HCRI70NTAL VELUCITV -IEl_3 U 
IM^LICI^ REALA8I A-h1o-Z) ” 

real * 6 PSI(Nl.NJ)tZ£rA{NI.NJ).UlNI.NJ).UP(NJI.AO(NOI.B(NJ)*ZCNJ). 
1 0(NJ) 

XDMM0N/CJ2DnLTA_X.J3XM2.JEj_TAZ,DZM2_,DZSQ.DZD6, j.Nl -U. NJMI * 1EXP2. _ 

I IS.JS.IE.J£,ISj».lEW,SP3Pl,JMAX,IS3.1=d 

UPI J I = ADI 1 ) 

DI I ) = 0. 

DO 3 I = IS.JE „ _ . 

aOUNOARY CONr^ITIO^S ON SECOND DERIVATIVE OF PSl 
03 = ZETAII,!) * DZD6 
DT = ZETA(I,NJ) * DZD6 

_CALCULATE..B MATRIX ElEMEO^XS 

DO I J = 2,NJ 

D(J) = (PSKI.JI - P5l(t.J-I)l X OELTAZ 

1 BiJ-l) = 3. « (DtJ) + 0(J-in 

J3I 1 >__=_BIJ J33 _ 

BINJ) = 3. * tOlNJ) 4- DT) 

TRIDIAGPNAL SPLINC matrix INVERSION SCHEME FUR J 

Z( 1 ) = BID 

J3U-X^-_= 2.NJ 

UL = -1 . / U=I J-I ) 

JPIJ) = AD(J) UL 

2 2(J) = B(J) + UL ♦ ZIJ-l) 

- 0 I I .N J)_ = - _Z(NJ)_^_UP-I-NJ) _ _ 

00 3 J = 2,NJ 
K = NJ - J ♦ 1 

3 JII.K) = - (ZIK) + Ull.X+m / UP(K) 

_ RETURN 

END 

tSU3RUUriNF, CALCLU PSI ,V.',FrR,FTl.AK,NI.NJ.NIM2) 

THIS ROUTINE CALCULATfS THE VERTICAL VELOCITY Ft£^J A 

IMPLICIT REAL *31 A-H.O-Z ) 

REAL* 8 PSIINI .NJ).v.(N1.KJ),FTR(NIM2).FTI(N1M2}.A<INIM2) 

C CMMON/C IZDEL TAX, 0XM2 , JE L T A Z , DZ M 2 , OZ S G . DZ Do , G . Wl “1 1 , NJM 1 , I EXP2 , 

1 I S. JS. I Et JE . IS .V , I Lfc , MPbPI , UMAX, 133. led 

C CALCULATE FOURIER TRANSFORM UF PSI 

DO, 3 . J = . JS. J E 

DO 1 I = I,NIM2 
F TR( I ) = PSII I . J) 

1 FT I I I ) = 0. 

CALL FASTIFTR ,FTI_,NIM2, !rxP2,-l_) _ _ 

C calculate FOURIER COEFFICIENTS 0<^ H 

00 2 I = I . N1 M2 
T = FTR ( I ) 

,..FTRII.)_=_- XIJ)_*_FIill). 
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2 FT I { I » = AK( I 1 • T 

CALCULATE W 3T TAKINi INVERSE TRANSFORM 
CALL FASTIFTR .FTI IEXP2, 1 ) 

W 1 N U±X * J L_= _fLI R { 1 J 

W(N1 . J) = FTRI2) 

DO 3 1 = IStNIM2 

3 W( I • J) = FTR( I) 

RETURN^ 

'ENP 

fSJ'^RQUTlNE FA| (PSI .ZEIA.F IRtFT I tPSl DFT . AD. BK.Si.ZR.Zl.JP.NI.NJ. 

THIS routine solves the POISSON EOUATtON F,^^ TrIE S r RE A MF INC T 1 PM 
US I NS A FOURIER TRANSt-ORM IN X. »^ARCH1NG bULJllClV IN ^ TFCHrJIO.IF 

1 N1M2.NJM2) 

_ IMPLICIT REALA3I AtH.O-Z.) 

REALAS PSI(Nl.NJ>,ZETA(Ni,NJ).FTRlNl.NJ),FTIlNl,NJi. 
i PSIBFT(N[M2.2).AP(NIM2|,RR(NJ).ai(NJI.ZR(«j|,Zi(NJ)«UP(NJ) 

C OMMON/C 1 /DEL TAX, OX M 2 , DE LT A 2 , DZ Mr . OZS Q , 0Z06 , S , .xl 4 I , NJHl . I EXP 2. 

1 I S. JS. I Ej JE , 1 S* J I EW.NPBPI , UMAX, I SBjlI Ed 

CALCULATE FOURIER TRANSFuRM OF ZETA 
OO 2 J = 2.NJM1 
DO 1 I = 1.NIM2 

FTRII.J ) = Z E TA( I . Jl 

1 F T K i . j I = 0 . 

2 CALL FASTI FTR < I , J ) ,FI I ( 1 , J ) ,NIM2. IEXP2, 1 > 

CALCULATE B MATRIX ELEMENTS 

_ _ ao 5 J_=_I.NIM2 

aO 3 J = 1.NJM2 

BR(J) = - ozsa * FTRII.J+I) 

3 RI(J) = - OZSQ * FTMI.JFI) 

BH( N JM? ) _=_BRJNJ IT ) _t_PS IBFTU ) 

3IINJM2) = BIINJ4?) + PSIdfT(I,2) 

CALCULATE FOURIER COEFFICIENTS OF P SI BY TRiJiAOONAu MATRIX 
INVERSION SCHEME 

UP! 1 J_=_ADITJ -- 

ZRI 1 ) = BR( I I 
ZI ( I » = BI ( I I 
DO 4 J = 2.NJM2 

UL /_ypLJrXJ - . _ 

UP( J 1 = API I 1 4- UL 

ZRIJI = QRIJI - UL * ZRIJ-n 

4 ZIlJ) = BIIJI “ UL A Zl(J-l) 

_ F TR < I ,N JMl ) _=__Z_R( N JM2 ) y _UP <NJ 42J 

FTIII.NJMII = ZIINJM2) / UPINJM2) 

DO 5 J = 2.NJM2 
K = NJM2 - J A 1 

FTRII.K + II = (ZRIK) A FT4(I,KA2n / UP(K) 

5 FTIII.KAl) = IZIIK) A FTI(I.KA2|) / UPIKI 
CALCULATE PSI BY TAKING INVERSE TRANSFORM 

._..DO_£> J = 2.NJM1 

6 CALL FASTIFTR (1 , J ) .FT II 1 , JI ,NI M2, IEXP2.-1 » 

RETURN 

JS UJROUT I NE F A STk X ,Y ,N. M. INV I 

THIS ROUTINE CALCULATES THE FAST FOUPIFR TRAXSEURM n- TH^ rnsiJtFX 
ARRAY (X.Y) IN PLACE, A'HrRE TLE I FNCHH n= X AN) F 1 ; A PIxFP nc 
J »0 

IM=>LICn REALA8I A-H,0-ZJ 
REALA8 XINI.YINl 

IMAX. =.N 

PIM2 = 6.2331 853071 BDAOO 
DO 2 L = 1 ,M 

JOELT = IMAX 

IMAX = IMAX / 2 

FJ = JDELT 

ASG = INV A PI M2 / FJ 
C - DCOS(ARG) 

S_= aSINIARG) _ _ 

U = I • 

V = o. 

DO 2 I = t,IMAX 

aO..I..J = I. N. JOE LT __ 

K = J A IMAX 
XJ = XlJl A XIKI 
YJ = Y(J> A YIKl 

XK_= X(J)_ji XJK) 

YK=Y(J)-Y(k) 

X|KI = U A XK - V A Y< 

VIK) =UAYKAVAXK 

XU) = X J _ 
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'I rcj) » yj 

T=C*U-S»V 

v = c«v-«-s*u 

2 U = T 

J = 1 

NT = N / 2 
I MAX = N - I 

30 5 1 = I.IM AX ■ ■ 

IFll. GB.J) GO TO 3 
T - X(J> 

XtJ) = X(I) 

-XiXJ—^.i 

T = Y(J> 
r ( Jl = Y( I » 

Y ( r ) = T 

3.K..= _NT 

4 IF(K.GE.J) GO TO 5 

J = J - K 

K = K / 2 

GO TO A 

5 J = J ♦ K 

IF(INV.EO.l) RETURN 
FN = N 
D06I=1,N 
X( 1 I = X( I > / FN 

6 Yt I I = Y<I> / FN 

RETOaJN ■■ - 

FNO 

iSURRnUTlNr GFk X.DET.Nl ,N2. I ACC ) 

THIS ROUTINE CALC'JLATES THE Ot T F g) M 1 N'A NT OF THE >1 A H 1 X X. TN oiflCF 
JSING GAUSSIAN ELlMlsAliPN w 1 T -l .^AkTIAl.. P 1 J T I 

IMPLICIT REAL*3{ A-H.3-Z) 

..REAL»^_XiNljbLl_l 

N2MI = N2 - I 
I ACC = 0 
OET = 

. _Q0 7.J^.l.N2MJ 

FOR EACH ROW, DETERMINE LARGEST MATRIX ELEMENT 
XMAX = 0. 

K = 0 

DD.l l_=^J..N2 _ 

Y = DABSIXI 1, JJ I 
I F I Y.LE. XMAX) GO TO I 
XMAX = V 

’ r CONTI NUE 

IF(K.NE.O) GO TO 2 
OET = 0. 

RETUR N 

2 IFtK.EO.I) GO TO A 
INTERCHANGE COLUMNS IF NECESSARY 
DO 3 J = I . N2 

T_=. XI J»I) 

XI J. I ) = XIJ.K) 

3 X{ J.K) = T 
OET = - DET 

IPi_J= !_>_! ... _ 

MULTIPLY ROW CONTAINING PIVOTAL ELEMENT BY APPRU^WIATt CONSTANTS. 
ANO ADD TO PDWS BENEATH, SO AS TC ZERO OUT AlL £_cM£NTS i3 = NEATH 
PIVOTAL element 

.00 6_.J..=. IPI.N2 . . 

I F( X ( J. I ) .EQ. 0. ) GO TO 6 
3 = XIJ.I) / XI I, I) 

00 5 K = IPI.N2 

. . 5 XI J.Xl _=_JU J) _l!_XCI.XI. 

6 CONTINUE 

7 CONTINUE 

calculate determinant FROM DIAGONAL ELEMENTS 

DO a I = 1 . N2 

I F ID AB5I OET » . LT. 1 .0^50 > GO TO 8 
DET = DET ♦ 1 .0-50 

1 ACC = I ACC F I 

. 8 DET ..= _0 El I 

RETURN 

END 

ts UBRQUT I nE GR at )C) ( X ,A ,LI N.E3 A.Ll Nf 3. NG ,NLE V.Ni . N J .N JSMl ) 

THIS ROilTINF CREATE^ i^ SHAOLD. LINE PRINTER GRA-JiUCAl. PtSPl.AY 
OF THE FIElD X 1 I H I \ T H i-. 0 izF 1 ME LIMITS i;F l\T -t.-.T 
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implicit REAL*8( A-H.3-ZI 

REAL«a X<NtiNJ).A(M.NJI.L INE4BI I 4 },L INES8( 14 > ,LIN^6y( 14) . 

1 OLn( 13 2) .NEm( I 32) . VI < 1 3 ) . V2( 13) 

__REAL44 X1NE14I 33) »LI_NE24t 33) .LINE 34 1 33 .NJ SMI) .NA 4 £Ll2j6J 

LOG) cal*) LINE2C 1 32 ) . L INE3( 132.NJSMI ) . SVH (33) 
equivalence ( LI NE2< 1 ) .LI NE24( 1 ) ) 

C OMMON/C) /del TAX«OX42.DELTAZ.DZM2.DZSO.DZ06.G.NI4I.'4J41.IEX0 2, 

» IS. Js.iEj jL=,is>:.iex.N?»uPi, UMAX, isij. lea. 

Ci3MMDN/C2/LINE4e»LINE58.LI Nhfc S , L I NE 1 * . L INE2 4 . N AMe, 

1 MT.IN.JN.INMl.NIS.NISMI.IL.SYW 

Z DETERMINE MINIMUM AND MAXIMUM VALUES IN THE FIElJ. A^J STORE FIELD 

r URSlDE_DOI»N_JLbmi_A 

XMIN = X(IS.JS) 

XMAX » XIIS.JS) 

OD 5 J « JS.JE 

K = JE - J ♦ I 

30 5 I = IS. IE 

V = X( I , J) 

A( 1- I S*I ,K) = V 

I F( Y.LT. XMINI XM l K = V 

5 IF(Y.GT.XMAX) XMAX = Y 

C CALCULATE INCREMENT BETlEEN SUCCESSIVE LEVEei UF SHADING 

D = (XMAX - XMINI / NLEV 

_I.F3D.GT^) .D-L2) .GD. 1J_XD - _ — . 

WRlTE<6,200t) <NAME(I,Na).I=1.2).MT.XHlN 
203 1 FDRMATI IHl . 2A 4, • FIE_D. TIME STEP f.I4.*. IS VIRTJAi_i_Y CONSTANT*. 

1 1 OX, "CLNST ANT =*.IPD10.2) 

RETURN 

10 OL = DLUGIO(D) 

IDL = DL 

IFOL.LT.O.) IDl =.^10L - 1 

SIGPRT = 10. »♦ ( DL-l OL ) 

IF(S IGPRT.GT.2. ) GO TO 15 
SIGPRT = 2. 

GO TO 25 

I A._I F ( S I GO^T . G T.S. ) GO TO 20 

SIGPRT = S. 

GO TO 25 
20 SIGPRT •= 10. 

_2.S X INC SIGPRT _*_1J),__**__IDJ ..... 

C calculate MINIMUM AND MAXIMUM VALUES TO 0E PwUTTHJ 

DV = (NJSMI ♦ I ) * XINC 

IMIN = XMIN / XiNC 

I F( XM IN.LT.O. )_1M I N_= LMIN - 1 

XMIN = IMIN ♦ XINC 
IMAX = XMAX / XINC 
I F( XMAX.GT.O. ) IMAX = IMAX * 1 

XMAX =.1MAX_4 XINC .. . - 

C 4RITE HEADING 

ih RI TE (6 ,2002 ) (NAME(I.NG),l = 1.2i.MT,XMIN,X4AX.XiMC.lS,)3.LlNEl4 

2002 F DRMA TI IHl , 4X .2A4 , • .-lELD. TIME STEP =*,I4,*, =*,1P010.2, 

1 . «. MAX If)UM =_'.J.1.3D10.2..*J_JINT£R VAL_=* ,lPD10.2x_ 

2 •, BASE POINT = < • , I 2, • , • , 12. • ) •// lx. 3 JA4 ) 

DD 50 J - 1 , JN 

LINEAR INTERPOLATION TO DETERMINE THE VALUE DF Trie FIELD AT EVERY 

? R l.N I_.CH AR AC T E R 1 N. A L IN E - . 

DO 30 I = I.INMl 

DIFI = (A(I*1,J) - A(I.Jl) / NIS 
L = NIS ♦ (I- 1) +2 

NEWIL) = A( I, J) - XMIN 
DO 30 K = I .NISMl 
30 NEMX+L) = NE'*(KFL-I) *• DIFI 

N E W { I L)_= _A.C I N. J ) - X MIN 

I F( J. EO. 1 ) GO TO 40 

LINEAR INTERPOLATION TO DETERMINE VALUE OF FlEe) uN ^INES BETWEEN 
VERTICAL GRID POINTS 

DO 35.1. = .2 . J L 

DIFJ = (NE.>(I) - OLOID) / OV 

TEMP = OLO(l) / XINC *■ 2. 

D J 35 K = I .NJSMI 

r £MP_.= _JEMP_ J-_..D1EJ 

L = TEMP 

DETERMINE SYMDOLS ASSCCIATED WITH THE VALUES OF Trte FIELD AT EACH 
PRINT character. AND =»hint LINES 

35_ L INE3 I I . K ). = SY.M (LI . 

wRl IE(6,2003) LINE34 

2003 FORM AT ( 1 X.33A4 ) 

40 DO 45 I =2,IL 

_L. = .J>JEW( I) / XINC * 2 . 

45 L INE2 ( I ) = SYM(L ) 

WR(TE(6,2003) LINE24 
DO 50 I = 2.IL 

50_OLD(.LJ._= .NEW.LJJ 
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MRITE(6*2003) LINE14 

C CALCULATE THE VALUE 3F EACH LINE SEPARATING l_EVii_S OF SHADING 

E = OMAX 1 (DABS( XM IN) , OABS( XMAX> 1 

EL. =_PLCGiQIE) _ 

I EL = EL 

IFIEL.lt. 0.) lEL = tEL - 1 
V2<I» = XMIN / 10. *» ILL 

VINC = 2. * X1NC_ y 10. ♦♦ l Ei 

VI ( I ) = V2( I ) ♦ V INC 

Vine = 2 . ♦ V I NC 

DO 55 I =2.13 

_ VUI) = VlII-1 ) » VlNC 

55" V2< I ) = + VlNC 

C WHITE SHADING VALUE SCALE 

WR1T£(6.2 00 4) VI. LINE 46. LINE 53. lEL.LI NE66 . V2 

__ZQ0 4 _EaRHATI DHO. 15X. 1 3E6.2 /JX . • SCALE 1 TO BE ’.1 AA 6/ J X._" jAJliJI PLJ ED 

1 lAAS/lX.'BY 134«« . I3.2X. 14A8/I 2X.1 3F6.2) 


RETURN 

J..ND 
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GRA'^H^ ( LINF3 4.LINF3.NJS MH 


THIS ROUTINE INITIALIZES ANNJTATION FOR GPA^HICA, OUTPUT 

__R£ALT.fl J_iNEABX-L4) «1 IVLSBiJ » J j 1_IN£C3TJ_A J .3LAN-ifl 

REAL *4 LINE 14(33). LINE241 33).L1NF34(3 3,NJSM1) .NA4E1(2.6). 

1 NAME2(2.6) 

logical*! LINE1(132>.LINE2(132>«LINF.3(132.NJSMi).L!NE4(112). 

a LJNE51 1 12 ).»_JNE5I 112) . SYH 1 1 5 3 ) . Sr.421 :>3 ) > O J1 .DASH. BLANK 

E GUI VALENCE (LINE1(1).lINE14(1)).(LINE2(1)._1NE24(1)). 

1 (L1NE4(|),LINE48(|)).(L1NES<1).L(NE33(1)>. 

2 ( LINE6( 1 ) .L I 'JE66{ 1 ) ) 

C OM MON/Tl /-.OFL TAX. DX M2 . OELT AZ. D2 M2 > DZS Q. DZDS, •la.j'tJ » J.aEX.PZ. 

1 IS.JS.IE.JE.ISa.IEW.NPBPI.JMAX.ISB.IEB 

C OHMON/C2/L INE4«,LlNt58._ INC6P. L INE 14 .LINC24.NAME2, 

1 MT.IN.JN.INMl.NIS.NlSMl, IL. SVM2 


DATA 

SYMl/* • 

. • a • , • 

• . *B* . • 

■ , *c • . * 

• • *0* .* •»•£*.* •.*£* 

.* • 

.* C* 

1 

• • 

,*H* ,» 

• .* I • •• 

• . • J* . • 

• . *K* »• • . *w* •• • . *M* 

.* * 

. * N* 

2 

• • 

,»o* .* 

• . *p • . • 

• . *Q* .* 

• .*R* .» • .* S* .* * .* T* 

.* • 

.* U* 


• • 

. • 1/ 1 . • 


•X • 




data 

NAMEIZ* 

r:* 

PSI • . ' 

• . • RHO* »• •.•ZETA*. 



1 

■ 

• • • 

u* . • 

• * • 

W*.* LOG* .‘(RIJ */ 



D ATA 

blanks/* 


• / . BL ANK4/ • 

•/.BLANK/* */.0ASr1/*| 

*/. 



_ a .DOT>’,»_.EZ 

CALCULATE NUMBER OF PRIM CHARACTERS FROM ONE H3<lZClNTAL GRID 
POINT TO THE NEXT 
MS = 129 / I NMl 

_ NFB = J.2 » N I S * INMI i / 4 » 1 ..... 

NI SMI = NI S - 1 
IL = NIS * INMI 4- 2 
INITIALIZE TITLES AND SYMBOLS 

DO 1 j = i^e 

1 NAME2(I.J) = NAMEKI.J) 

00 2 I = 1.53 

_ 2. SYM2ia)_=_ SY'41LIJ 

initialize arrays used in PHINT LINES 
DO 3 I = 1.33 

3 L1NE14(I) = BLANK4 

_ DO 4 i__= 1 . IN 

4 L INE I (N I S*( I- I ) *2 ) = DOT 

00 5 1 = NFB. 33 

5 LINE24(I) = BL4NK4 

L1NE2 ( 1 J_=.DQI ... 

LINE2(ILFI) = DOT 
OO 6 I = I.NJ5MI 
L INE 3( 1 . 1 ) = BLANK 

DJ 6 D.--NFJ3.33 ._ 

6LINE34{J.I) = BLANK4 
00 7 I = 1.14 

LINE481I) = blanks 

L I NE 5«( I ) _= .BLANKS 

7 LINE 68(1) = BLANKS 
OO 8 I = 1.13 
L INE4 ( 8* I + 1 ) = DA SH 

a _ 1NE6 (8*1-3 )__=_D.ASH .. .... _ _. . 

03 9 1 = 1.26 
9 L INE5(4*l4-2) = SVM1(2*I> 

DO 10 I = 1.52 

_ 10 LINEd(2.*1f3X_= dash 

RE TURN 

I.swoRMjT K.r pr5'~*Er.} ( irrpE.ia.ja.Ni.Rj . npbpz ) 

THIS ROUTINE OFTERMINFS THE TY^F OF EACH GRIP J,y 1 N ( IN T.HC FIELD 
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implicit REAL*8< A-HtO-Z> 

.JNTEGER*4 i I \T P£1 N I_.NJ J 3 INPJ3R21 ^ JBI NPBP2 J 

CO'4MUN/Ci/’DELTAX,CXM2.DELTAZ.C>ZM2.DZSQ>DZD6»tf.Ni<«l.NJHl.l EXP 2, 

1 IS«JS.IE.Jc«lS».IEh.NP3P|,JMAX.ISB.lE3 

C initialize all internal POINTSf AND UPPER ANJ 3uUNDARIES 

Da_2 X s . X*J^ . _ 

DO I J = 2*NJMI 

1 ITYPEII*J» - 0 
I TYPEI I .NJ) = 1 

2 I TYPE (1*1) = 2 

determine type of all points on the SURFACE UF THE BARRIER. 
PROCEEDING FROM LEFT TO RIGHT 

ISO e__IB{2_) 

lEB = IBINPBPII 
1 S( I I = ISO 
J3( I ) = 1 

_ _l a< NPBP2 l_=_J EB 

J0<NPBP2) = I 
I TYPEtl SB. I) = 5 
I TYPEI lEB.l ) = 6 

DO 9" I = 2,NPB(^ 

K = IBI 1 1 
L = JBI I I 

I FfL.GT. JMAXl JMAX = I 

IF(JBlI-l)-L) 3.4.8 

3 ITYPEIK.L) = 3 
IFIL.EQ. JBI I4 t n ITY’EIK.L) = 7 

GO T0_9 

4 I F I L- JBI IFl » ) 5 ,6 ,7 

5 I TYPE! K.L) = 5 
GO TO 9 

. 6. JTVPEIk_.JL1_=._2 

GO TO 9 

7 ITYPEIK.L! = B 
GO TO 9 

a ITYPEIK.LI = 4 

1 FlL .Ea.JBI IF I > ) ITYPEIK.LI = 6 
9 C ONT 1 NUE 

JMAXMI = JMAX - 1 

DEFINE_ALL PO I NT S J N_S_I D_E._T HE B A_RR I EP TO _BE TYPE 9 

bo 10 J = i. JMAXMI 
K = 0 

DO 10 I = ISe.IEB 

_L = I TYPE I I, J) 

1 FIL .LT.'3.AND.K.FQ. 1 ) ITYPEII.J) = 9 
lFIL.EU.3.0R.L.Ea.5> K = 1 
10 I F ( L . EQ. 4.0R. L. EO .6) K = O 

RETURN 

', £NP 

15UL3ROUTINF S^ ( X . P SI AUX . 1 T Y PE . BET A . I 0 . JO .X 0 . N I . InJ . N I MZ . N J M 2. N IP A . 

I J 'lA XMl , NP3 . NPBP2 > 

T~Ri~S~RPUTINF 5UPE nPQ'SCS~^'HE SDLUTro>^OF~~T HE' ~^i'5~S'uV' U J A f 10N~wTfH 
THE EFFECT ON PSl OV THR VU T. T 1 C I I Y i-Rj-1 EACri PUINT ON 

THE 3UI<r ACE OF the OAkRIEI* 


IMPLICIT REAL*3( A-H.O-Z) 

REAL *8 X{NI.NJ).PSIAJXCNr?A,NjM2,JMAXMl).D£TA(NPJ.NP3l 
I.NTEGERF4 ITYPE(NI.NJ),IBI NPBP2 ) . JB I NPHP2 1 

C py.MON,^Cl /DELT AX, Dx:A2 . OE LT A Z , DZ M? ,DZS 0 .OZDS.i^Nl'M.NJ'U .1 EXP2 . 
i 1 S. js . I E. JE . IS A . I CX .NPOPl . JM AX, ISd. lEP 

03 2 K = l.NPB 
ALPHA = 0. 

D3 1 _L =_1.NPB 

C CALCULATE ALPHA FROM BETA 

1 alpha = ALPHA + 3FTAIK.L) ♦ 1X0 - XI I O ILf 1 ) , J d I - F 1 ) } | 

L = IBINP3PI1 - XeiKFl) 

M = JRIKFII - I 
DO 2 I = I . NI M2 
03 2 J = 1.NJM2 

2 X I 1 , J FI |_=__XI J . J* 1>__F_ Ai__PHA _.*_,PSlAUXt L*I , J.MJ 

C DEFINE PSI AT SlOE BOUNDARIES 

03 3 J = 2.NJ-4| 

X(NI’-II.J) = XII, JI 

3_.X ( NI , J) = X ( 2 , J) 

c redefine uoundary condition exactly cn and inside barrier' 

OO 4 I = ISO, lEB 
DO 4 J = 2. JMAX 

4 I FI I T YREtl. J 1_.NE._0I_X( J . J 1_=__X0 

RETURN 

END 


89 



Listing of data cards 


65 

24 

ao 

to 

0 

0 

6. 25^0402 6.250402 

t.0004Ol 

12 

1 

2 

1 

l 2 
65 

3 

24 

1 2 
2 

4 

12 5 13 5 14 

5 15 5 15 4 

40 
60 
80 
3 0 
80 

05 -P 

00 

25 . 

25.. 

.25 

..9 




80 

80 

2 

2.sor>4-oi 

1 

0 . 0004^0 0 2 . 500 402 


3 15 


2 


Beginning of output 


SUM«4A^f OF INPUT DATA FUR NON- 1 T£R AT 1 V£ • SOLID tiARRILK PKGORAN 

Nl s NJ s 24» NT K aO« NPO ■ 10. iiJT > U« ISAVE >0* DX 6.26J402* OZ « 6.2SD4^02» DT = UOOO+dl 

'^COOR"Dl\ATES"'OF~aAHRlER“ie' 'l 12. 21 i I273t T12.VI " ~l 1 2 . t» * ' i 1 3. 5i ^ll4'.5"|i UVtSJ U5.A> U5.3) ”ll5ia> 

IS ■ i. JS = 1 . IN = 6S« JN > 24. ^JSMi * 3 

NPSILV » ’"40 ,”^ NRMOLV^ 40 i NZIALV"”^ 2S"i NULV""*' 2S. ’ NWLV "zs", NRlLV 9 

IPSIGR >60 50 000000000000000000 

IRH3GR ’ 80 0 ~ 0 0 "" 0 0 O'^^O 0 0 0 0 " 0 0 O”’ 0 0 ~ 0 ” O' ' 0 

IZT4SR > 80 0000000000000000000 

iuGR >"”80 "“'O 0 “ "O O 0 ' 0 O ' 0 ” O ‘ ' 0 0 0 ^0 0 ‘ O 0 0 0 ~ ~o ^ ■ 

l«GR > 80 0000000000000000000 

IRIGR « 60 ’ 0 0 ~^0 0~^0 0 ' "O ’ "O .0 0 0 " 0 0 ' 0 ^ O" 'O 0 ^O' 0 ‘ ' 

CONSTANT V£-QC1TY CASE. U1 ■ 2.50040 1 

CQNST4NT LAPSE RATE ‘case* ~ gamma' » 0.0 TO » 2. SOD 402 
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k' 


' — PSl'|C|a_D. TIPE STEP ■ *0. MININUM • -S.kOO'pM. wVi'miM • 1 1 NfcHV*l. • 1.600*04. BASE POIPT • < I. I I 

.A*I***i***i«»A*A*AAAA*»*»***A*44*44»*AA**AAAi*4AA4AA»AAAAA/Ul/!AAAA*AA/L»4AAAAAAAAA4AM«AMAAAAAAA*4»****»*»*»**»»*‘*« 
AAAAA4AAAAAAAAAAAAA AAA AAA A A AAA AAA A 

BBBBBBBBBBBBBBBBBBr ®®®®‘’®®®®“““®®*“‘““®"““®“®®“®““®®®““®®Si2aBaSBlHBiliBBBBBBaBB8BB80SflB 

BBSBBHBBBB88BBB CCCCCCCCCCCCCCCCCCCCC BBaBBBSBBBBBaBBB8BBBBBB8BB88SB8 

AAnMAAAAHA CCCCCCCC^CCC CCCCCCCCCCCCCCC-- - - - BB00BI?ooOeBB9PtfWHPOO»»i^i. 

— ccccccSc CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC^ 

CCCCCCCCC OOOOOODODOOOO ■' **'■ ** 

cccc cccccc ddoddoododooddooddddddod 

.* CCCCCCCCCCCCC 0000^00^ DODOODDDO^ 

ccc-c:c DDoooSo°°° eeeeeIIeeeIeeeeeee ooooooDODOoooooODOoboooooooDoboBBoDODDODD cccccccccccccc::c:cc 

.X ......i........ oonr^nnonoDDOOD ccccccccccccccc 


B8QaaSB9Bda)39»B9 

CCCCCCCCCCCCCCZCCQCCCCCCGCCCCCCCCCCCCCCCCCC'e 

cccccccccccccc 

CCCCCCCCCCCC • 

0OO0ODOO00DDODDODO0D9ODOOD0OOOO00O3 30O •. CCCCCCCCCCCCCCCCCCCCCCCCC 

E€eEEEeEE££eEeCEE£ doooooooddooooddooooooooooodddoodododdooo 

eeeeeec eceeecec ooooddddoooooo 


CeEEECEE 


^ DDODDOO - 

. OBDOODD EE£EE£ EtECEEC _ 

ODOOOODOO EEEee ■ 'FFFFFFFFFFFFF CEEEEEEEeEEEEeCEEEEEEEEEeeEEEEEEE 

0000090 DOOOO EtEEC FFFFFFFF FFFFFFFF E CEEEEEEEEEEtEEEi ECEECEEeeEeEeEEEEEEEE 

00000300D EEEEEE FFFFF FFFFFF 

A.QQQQ _ _ E£E££e FFFFi^ GGGQCG _ -FFFFF 

’ EEEEEE FFFF GCCCCCGGGGGCCQ 

EEEEEEEE FFFF CGGCGC GGG6GG 

EEEEEEEEEEC FFFF GCGO QGGC 

>EeEEEl5Ee...-_-^ FFFFFF 
eeEEg FFFFFF 

FFFFFFF 


DOOOOOOODO 

OODOOOOODO 

00000000000 

EEEE EEEEE EEEEE EE OODOOODOOODODOODOOODO 

EEEEEEEEE OOOOODOOODDOOO JOO * 

EEEEEEEEE DDODOO 

^pp^PfrpFFFFFFFFFFFFFFFF FFFFFFF* FFFFF EgEE EEEEE 

ffffffffffffffffffffffffffffffffffff eeeEEEEEEEE 

GGGG HHHHMHHM GCGGG F.FFFFFFFFF ... EEEEEEE EEEEEEEEEEEEEE f. 

GGGG HmHMMMHHHHMHMH '* G6CGGG ” FFFFFFFF EE£ E EEEEEEEEEEEEE 

GGCC HHHHH HHHHH GGGGGGGGCGGCGCGGGCCGGGGCCCCGGGSGC FFFFFFF 

AFFFFFFFF GGGG HHHHH HHHH GGGGGGGG GGG GOGC GGGGGG GGGGGGGC GGCGGGG FFFFFFFFF 

■ FFFFF'ffffeFF.. . GCCGO hhh* 1_ .UlMItl _ HHMhH _GGGGG - * ^^^^ffffffffJfffffffffff*" 

fffffffff cgggg hhh 1 1 j 111 I It 1 1 1 1 HHMHH 

FF CGGGGGG HHMH UMI I III MHHHHH HHHHHHHHHMHHM GGGCGGGG FFFFFFFFFFFFF 

GGGCGGGG MHHH till III! MHHMHHHHHMMMHMMHKMMHMHHHHMHHMHHHMHH GGGGGCGG 

.ii" "5! 

ml JJJ KK . JJJJ laillllUIIlllJilllULUIlLLlIIIll.- hhmmhhhmhhhhmh 


HMHHHHHHHMrt 


I m 
mil 

It II 1 1 1 

iiiimiu.. 


mil j mu 


fill 


HKKKKKKH JJJ 
KKKKAK<AKKA JJJ 
KKIC KFKK JJJJ JJJJJJJJJJ 

KKK .AKA _JJJJJJJJJJJJJJJ JJJJ JJJJJJJJJJ JJJ 


JJJJ 
JJJJ 
JJJJJJ 

- JJJJJJJJX 

JJ JJJ JJ JJJ JJJJ JJJ JJJ JJJ JJJ JJ FAK 

J JJJJJJ JJJJJ JJJJJJ JJJJJJ KKKK 

KKKK 

• KKKK 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


mill 
mil I til 

miiirim 

_ _ Lniuimium * 

HKK jjjjjj JJJJ j JJJJ jjjjj JJJJJJ jjjjjjjjj iiimiiiiimiMiiiMi 

AKA JJJJJJJJJ JJJJJJJJJJJ IllIlllIIII 

KKKK JJJJJJ JJJJJJJ 

KKKAAKKAKKKKK JJJJJJJ JJJJJ JJJJJJJJJ - . - -A 

_ ^ i.LL KKKKKK KKKKKKKKKKKKKAKKKKKKKKKK J J J J J J J J J J J J J J J J J J J J J J 

LUU MMMHM LLW KKKKKKKKKKKKKKKK KKKKKKKKKKKKKK JJJJJJJ 

MHHAMMFM U-L kkkkkkkkk KKKKKKKKKKKKAKKAKK 

LLL MHHMA MFMM LLL 

LL NMMM MMH LLL LLLLLLLLLLLLL 


KK K LLLLL 

AKA LLLLLLLLL 
LLLLL LLLLU 


• LL . LL . 

LLLLULLLLLLLLLLLLLLL i-LL 

LLLLLLLLLLLLUULL LLL HHM 

LLLLl-LLLL LLL **•* 

LLLLL CLLL km NI 
MMHMMMHmMMMMMMMMHMMM LLLLLLLL MM NI 

FFMMMHMMM LLL **AM NN 

MMMHHH MH NN 

•NNNNNNNNNNNNNNNNNN hnhmm HMM nn 

NNNSN'4NNNNNNNNfF6NNNNNNH M>4F HMM NN 

NNNNNN AHMM MMMM NN 


kkkkkkkkkkkkkkkkKkkkkkkk<k<k< • 

kkkkkkkkkkkkkkkkkkkkkk 

LLLLL LLLLLLLLLLLLLLLLLULLLLLL KKK 

LL MM NN MM LLLLLLLLLLLLLLLLL LLLLLLLLLLLLLLLLLLU-LLLLLL 

LL MM NNNNHN MMM llllllllll llllllllllullllllllllclllllllllllllllll * 

NNNNNNNN MMM LLL 

NNN Nfj4 MHM HMMMMMMMMMM 

NNN NN MMMM MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM 


LLLLLLLLLLLLLLLLLLLLLLLU 


NNNN 

«OOaOQOOQQQQOQaOQQQQOO N( 

OaOQOOOOO OOOOOO ' 

3 OOOO 

ppppp;>ppppppp 


HHMM MMMM 
FMHH MMHH ^ 

M MMHWMMM 


MHMMMMMMMMMMMM MMMMMMMMMM MMMMM MMHMMMMM M MMMH MMMMMMMMMN MMMMM MMMMMMMM MM MMM • 
00 NN MMMMMMHMMMMMM MHMMMMHMHMMMMMMMMMMMMMMNHMHMMHMMMMHM 

QOOOO NNN MMMMMMM 

OOOOQOO *^NNM .jnnnnnnnnnnnnnnnmnnnmnnnnnnnnnnnnnnnnnnnn nnnnnnn* 

00 OOO NNNN NNNMNf'WNNNNNNNNNNMHNNNNNNNNNNNNNNNNNNNNNNNMNNNNNNNNNHNN'JNNNN 

O 000 NNNNNNNNNNNNNNNNN NNNNNNNN ^F4NNNNNNNNNNNNNMNNNNNNNN'i'«NN 

10 OOO NNNNNNNNSNNN 

)0 FP OOO NNNNNNNN ___ 

I ppPP 000 OOOOOOOOOOOOOOOOOOOOOOOOOO OOOOOOOOOOOO 

) PPFPPP OOO □OQOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOODOQOO 

) ppppppp OOOO OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOaOOOOOOOOOOOOOOO 0000000 

pppappppppppppppppppppp OQ MMMMMM . N 0 pppPPPP 0 OOO OOOOOOOOOOO . .0 00 000 000000 O • 

pppppaopp PPP o N MMMM N O PPPPOPPP OOOOOQOQ003 

PPPPP FF oo N MM NN OO PPPPPPPP OOOOOO 

‘ QOQOOQOOQOOQOQQ P O NN NN OO pppppppp pppppppppp pppppppppppppppp 

uuuwu Jwuuuuuuu u n ppppppppp pppppppppppppppppPPPPPPPPPPPPPPPPPPPPPPPPPPPPC»P'*P*>PPP*»Fo • 

ppppppp pp P pp P PPP p PPP PPPP ppppppp PPP PPPPPPPt^PPPFP ppppppp PPP PPPPPPP 

pp PPP pppppppp ppppppppp 3 PPP PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP 

pppppppppppppppppppp pppppppppppppppppppppp OOOO 

pppppppppppppppp oqooooqo • 

pppp ppppp ^ OOQQOaOOOQOd 

ppp OOQOQOOOOOQOOaQOQ 

PPP QaoOOOOQQOOOOUOOOQOOOOOOOOO 

pp _ QQQQQQQOOOQ QQOQOOOQQOOOOOOOOQOOQOOOOOOQOOQOO • 

pp ' oQQaoQaoQQOQQooaaoiKjQOOOOoaoaooooQaQaooooooQoooQQoaoQoooo 
ppp QOQOQoaaQQaoooQoaoQoaaQQOoouQaooaoiiooooooQQoooooQ 
pppp OOOQOOQOOOaaOOQQOOQOOaOOQQOOQQOQOQQQOa 
pppp OOQQ QOOOQOQOOOOOOOQ ...... . 


OQO 


QQQOOQOOOQOO 

OOOOQOQO O OO D NN NN □□ 

OOOQOQQQQ RPPPMORRR P O N^J4NNNN OOO 

QQQQU RRRRPRRP RpQp O NNNNN OOO 

RARRRRHRRR SSSSSSS POO NNN 00 

SSSSSSS P O OQO 

SSSSSSSROP 00 OOOO 

SSSS59SR0 □□ OOOOO 

SSSS^SSRQ P oooooooooooo 

SSSSSSSRQ P OCOOQOOOO 
SSSSSSSR PP OOOOOO 

sssKsssn 

SSSSSSSR 


RRRPRRRRR 
R RRRRRRR 
R 4 QRRR 


• RRRP9RRR 
PRRRR9RRR 
RRRRR^ 
RRRRP 
*RRR 


sss 

SSS5 
SSSSSSS 
SSSSSSS 

ass sss sss SSSSSSS* 
sss sss SSSSSSS sssss 

SSSSSSSSSSSSSSSSSSS . 
SSSSSSSSSSSSSSSSSSS RR 


PPPPFPPPPPPPP OOO 
I ppppppppp OOOO 

a OOOO 

OQO OOQQO 

OQQQOQQQ 


RPR 

RRRRRR 

RARRRR 9 R 49 R. 

RrRARRRRRRRRRRRRRRRR 

RRRRRRRRRRRRRRRRR RRRRRRRRRRRRRRRRR 
PPPRR rrRRRRRRPRRRRRRRRRRPRRRRRHRRRRRRRRRRRRRR^ 

RR«RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRnRn«««»**'»R«'»««»^""'’*’‘’ - -• 

RRRRRRRRRRRRRARRRRRRRRRRRRRRRRRRRRRRRRRRRRRR 


SSSSSSSSSSSSSSSSSSS RRRRR RRRRRRRRRR 

* S SSS 55 SSS sis Is sss Islll Ills SSSSSSSSSSSSSSSSSSS sssss sss 953 SSS 1 S SSSSS sss St S 4 SS 9 S sssss sssss sss sssss S 3 3 SSSSS S 3 sssss sssss sssss SSSSSSS# 


SCALE I TO se 
MULTtPutEO 

OV lOP* s 


lo| |e< 

-2«se 


•60 
Im , 


2*20 

at IhI 


-l#«9 
111 \j{ 


k1 

•ilso 


c; 


D«60 

IpT 


|n| |n| 1 ol |pT " jal iRl l^srUr^' 

;2o * " “ ■“ " 


Tul Ivl |»| 


.00 

Ixl 
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PSI FIELD. TIME STEP •’~so7 MINIwipi • -J.E«D*»5. MXIHUH - l.'aOMO*. INTERVAL ■ l.•00*a•• EASE POINT ■ I I. II 


cccccccccc 

CO 

ccccccccccc 

CCCCCCCCCCCCCCCCC OD©DOOO 

_»ccccc:ccccccc. . ododpoo 

CCCCCCC 000000 

oooooo 

DODOOOOO 
DDOOOODOOaOQp 


aOBOBOeeeOBftbOBBOaoeBid30‘3 9B33B 
BD8Bb6eb6BB60BPB0BB5993«. 


ODOODOD 

. DQOOO 

FEEeeCCEEEE OOOOOO 


aaaaaa A aaaaaaa aaaaaaa aaaa aaaaaaaaaaaa aaaaa^a aa aaaaaa aaaaaaaaaaaaaaaaAa 
AAAAAAAAAAAAAAAAA AAAAAAAAA 

BBBBBJaBBBBBliMBla aSBBBBBBaBBBa BBBBBBaaBaflBaBBBBBBBBBBaBBBBBBBBBBBBBBBBB 

BBBaBBBBBBBBa CCCCCCCCCCCCC^CCCC^ CCCCCCltCCCCCCCCCCCCCC 

" - ccccccccc“ — “''"^'^‘^iccccSccccccci 

oaDOOODOOODDDDDO CCCCCCCCCCCCCCC ojoOODOODO ‘^“‘^“‘^CCCCCcIcCCCCCC CCC 

OODOQDQDDOOOODPODODOP . ___>^_CCCC.CCCCCCCCCCCCC CCCCCCCCCCCCC •_ 

DODOOOOO OODOOOOOOOO ^ CCCCCCCCCCCCCCCCCCCCCCwCC 

ee^eceee eceeeec ooooooooooooooo oooooooddoo cccccccccccc 

trppee EEEEE OOODOOOOO EEEEEEEEEEEE DOODDOOOOO 

Ieee wwwuui/ww EEEEEEEEOsEeeeeeeeEee ddpoooodoo t. 

EEEEeIe^^^^” eEeeeEeeE~ ODOOOOOODODOOOOODOODOOOOOOOOO 

pppIp FFP?? FFFF liiEEEEEeEEElEll EEEEEEEEE ODDODDDDOOOOOOOOODaDTDDSD 

OOODSJJOO peIIP' fV/f CGG PFF IeEEEEEEEEE FF-FFFFFFFFF EEEEEEEEE OODDDODDDODDOODa 

- - EEEEEEEEEE FFFF OOOGG ““=^“55 '"fFFF Ff'fF 0 '"'^ff'fFF '“IeIIeIeEEEEEEEEEEEEEEEEEEEEE 

FFFF OGg“ HHWMMMM “® 00G FFFFFFFFFFFF aSGGSaaCGGCG = 

FFFFF GGG HHHpNNHHHMMM ^GG FF£FF£fF_ GGGM4 # .._SGCGGCOGO PFFFFFF .. ..EEEEEEEEEEEEc= 

— FFFFFF “GGG MMHH “ MHHH GGG '^"JSJIJfffF 

fffffffffffffff'''' cgg“ hhh” iiiiiii **^h “gcggg GGOGG mnmnJIJh»<hmm GGGCGGO FFFFFFFFFFFFFFFFFFFFFFFFFF 

CGg“° HHh”_ IliliiLill SmM ®:GGGGGGGGQS NNIlHNHmHHNHHmHH GGCGWG rFFFFFFFFFF.FFFFFFFFFFFF^ 

’ GGg 2 gg“‘^ "hnS"" III" "l MNMM*"* HmHHHH “gGGGGGG 

,{{* JJJJJJJ ‘11 HMH HHNN IIIIIIIIII HHMMHMH GGGGaaGGGGGGGaaSGGGGSSGSCCGG 

HhhHMMHHh”*^ l!l* jH"* HKKKKK *^JJ * !l ””””” III” JJJJJJJ XlllII MHMHMMHMHHHMMHH 

mhhhhhhhhhh^IS^JSSh”.”^. iIlL^jUL 

nil JJJ KKK KKK JJ XllllIIXlll JJJJJ JJJJJJJ ‘MU!! 

mil JJJ kkk rk JJ iiiiiiii 

Him 11 JJJ KKK LL4.L KK JJ I III 
iinimimuiJiilt jjjj_jckk 

imiiiiiiim jjjj KA llll lll aa jjj 

JJJJJJ KK U. Li. 


JJJJJJJJJJJJ 


KAK 
KKKKK LI 

KKKKK<KK<AAKKKKKKKKKK<KK<AKKA U 

•kkkkkkkkkkkkkkkkkkkkkkkkkkkk LL 
KKKK LLL 


KKR LL M LL KK JJJJJJJJ KKKK 
_JLK dJJJ 

L KK KKK 

LL KKK KKK 

MM L KKKK KKKKK 
NNNN MM LL KKKKKKK 
NNNNNNN M LL KKKKK 


HMHMHMHHHHHMHHMH 

JJJ JJJJJJ XJXlIttll 

JJJJ KKK JJJJJJ txxMiiiimmi 

_ JJJ KKKKKKKKKK JJJJJJJ. * ^* ********** H I I H ? 

JJJJ KKKKKKKKKKKKKKK JJJJJJJ 

JJJJ JJJJJJ KKKK KKKKKK JJJJJJJJJ 


m 1111 m I X m 


-LLLLL- 


MN 


MMM 


JJJJJ JJJJJJJ JJJ JJJ 

JJJJ^JJJJJJJ/JJJJJJJJ JJJJJJJ fc, 

KKKKKKK JJ JJJJJJJJJ J 

KKKKKKKKKK 


LLLLLLLLLL 
LLLLLLLLLLLLLLL 
LULL LLLLLL 

LLL LLLLLL 

LLL M4MM LLLLLLL 

. I 1 I i HI i_L LL LLl‘" m” NN^ ^NN ^MH LLL LLl'" MMMMMMMmSmmmMM ^’"^'lLL^LLLLLLLLLLLLLLLLL LLLLL L LLLLLLL LLLLL 

LL.LL.LLLLLLLLLLLLLL^^^^^^^^^LLL ^M ^NN NN^ ^ L^L^ MMMMMMM LLLLLLLLLLLLU-LLLLLLLLLLLLLLL .L 

^.LLLLLLLLLL MM NN 000 MN MM LLLLLLLLL MMM MMMMMMM 

MMMMMMMMHMMHMMM LLLLLL MM NN OOOO NN MM LLLLLLL MMM '^'^MMMMMMNMIiMMMMMMMM 

MM MMM**iMMMMHM M NN OOOOOO NN MM LLL MMM MNWHNN HMMMMMMMMMMMMMMMM 

MMMMMM MM N OOOQODO NN MMM MMM MNNNNNNNNNNN 

NNNMNMNNNNNNN MHMMM MMM N 000 OOO NN MMM MMMM MNNMN^«<NNNNN^f'^N 

NNNNNNNNNNNN>J<NNN MMMMMMMMHMM N 00 0000 NN MMMM MMMMM NNNN NNNNNNNN 

- 000 NN MMMMMMMMM NNN 

00 NN MMMMMM NHN 


NNNNNNN MMMMWMMMH N 0 

^ □ NNNNN MMMMMMfc* N 00 

OOCOOJOO JOOOOOO NNNN mMMMMM N O 

0003 OOOOOOOOOO NNN mmmmm n O 

OCOOOO NNN MHMMM NN O 
^ pppspppppp 0000 NN MMMHM NN O 

pppoppppppppppppp 0030 N MMMM N 0 

pppp PPPPPPPP OO N MMM 

ppppp 00 N MMM 
• OOQQO pppop 0 NN MM 

OOOOOO OQQUOQOQOOQ PP 0 N M N 

OOOOOOO QOOQOOQOOOOOOOOO PP O N N 

00 00 000009 PON N 

«0 QOQ 0 NN N 

p pppppHR P OOP 0 N NNN 

RMqPHPPRPRRRRRPRRPPPPR O NNNNN 
RRPRR ROPRRRRRRRRRRRR RRQP NNNN 

•RRRPOiapR SSSSSSS o nnnnn 


MMMMMMMMMHMMMMMMMMMMMMMMMMMMMMM • 

MMMMMMHMMMMMMMM n 

NNNNNNNMN 

NNNNNM><NN NNNNNNNMNNNNNNN 

NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN 


NHNNNNNNNNNNNNKMNNNNNNNNNNN 


.=pp_ 


N O 

N OOOO 

N oooooo 
N 00 00 


0300 

ooooooooo 

000000000333 

000000000000000000 • 

0000000003 OOOOOOOOOO OOOO ooooooooooooooo 

OOOOOOOOOOOOOOOOQOOOOOOOOOOOOO p»p»»p 

OOOOOOOOOO PPPPPPPPP 


SSSSSSS 
SSSSSSS 
SSSSSSS 
SSSSSSSSSSSSSSSSSS5 
SS5SSSSSS5S SSSSSSSROP 


SSSSSSSSSSSS 
SS5SSSSSSSSS 
SSSSSSSSSSSS 
5 SSSSSSSSSSSS 

ssssssssssssss 

SSSSSSSSSSSSSS 

SSSSSSSSSSSSSS 

SSSSSSSSSSSSSS 

SSSSSSSSSSSSSS 

SSSSSSSSSSSSSS 


SSSSSS5RQ p OO 


MMM NNN 000 

OO NN NNN OOOOOOOO 

00 NNN NNN OOOOOOOOOOOO 

QO nn NNNN Q 0333aoaooaooaOQ 

00 NNNNNNNNN 00000 
OO NNNNNNN OOOO 

OOO NNN OOOO 9PPPPPPPPP 

nno OOO _ ... . 

030 ■ OOOO "PPPPPP pppppppppppp 

U PPP OOOO 00300 pppappppppPPPPP PPPPPPPPPPPPPPPPPPPPPPPPPPPPPP 

N 0 PPP OOOOOOOOOO PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPRPPP 

NH 0 pppp OOOOOOOO PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP 

NN 0 pppp”’ oooooo PPPPPP pppj^pppwpp 

0 pppp pppppp 

00 ppppp ppppp 

o . pPPPPP . . ppppp .... .. . 

0 ppppp ppppp Q QOOOQQOQ 

OO p pppppp QOOOOQOOaOOOOQQOO<30QO>«^..<..-^ 

QQ PPPPPPPP rkrkAAnAAi^nonnaaooopaQoaaoaaQooaaooooooooQoooooo 

oo pppppp-. 

OOOO ppppp 

OOOOO ppppp 


PPPPpRp PP 


sssssssRO p OOOOOOOO PPPPPP oaoaa 

SSSSSSSPO PP OOOOO PPPPPPP OQQO 

SSSSSSSR O PR PPPPPPP OOO RRRPRR 

SSSSSS5R 0 PPPP PPPPPPPPP OOOOO RRRRRRp RPRRRRRR PR 

SSSSSSSR QO PPPPPPPPP QOOOQ 

SSSSSSS R 00 QOOQ 


QO 
900 

OOQQOOO* 
OOOOOOOOO 
009000930 
OOOOOOOOO 

OQQOOOOOOOQ • 

QQOOOOQQOaOOO R 

0 00 OOOO 00 00 OO OOOOO 0 00 OOOO 00 QOQQ OOOO OOOOOO OOOOOO RRR 

QOQ OO OOOO OOOO 0000009 OOOOOO 90000 00 OOOOOOQOOaoOO 
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LEGEND FOR FIGURES 


The figures show the line printer output at fixed times for selected 
fields in the (x,z) plane, that is, in a cross section normal to the axis of 
the disturbing obstacle. The top line contains the symbolic name of the 
scalar field, the time step number, the range of the values plotted, the 
interval from one band to the next, and the grid point (i,j) values for the 
lower left hand corner of the output field. Value intervals are represented 
by alternating bands of letters and clear space. The values of the equi-scalar 
contours constituted by the boundaries of these bands may be read from the 
scale at the bottom of the diagram. The scale is subject to change for each 
diagram, since it is automatically adjusted to the range of output values. 

The grid points at which the scalar field is calculated are delineated 
by dots around the edge of the diagram. Values between the grid points are 
determined by an interpolation routine. The barrier is indicated by a 
solid line at the bottom of each diagram, except for figures 13 and 14. For 
this case only, the lower boundary was specified by a nonrigid flow condition 
located between 11 and 20 grid points from the left side. 
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Figure 1, Streamfunction field at 1000 seconds for case 1 (constant 
velocity, constant stability). Ax = 1000 m, Az = 1000 m, 
horizontal extent displayed = 64 km, vertical extent 
displayed = 23 km. 
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Figure 3, 


Streamfunction field at 2000 seconds for case 1 , 
as in figure 1 . 
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Figure 5. 
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Vorticity field at 1000 seconds for case 1 


in figure 1 . 
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Vorticity 


field at 2000 seconds 


for case 1 . 


in figure 1. 
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Figure 8. Richardson number field at 2000 seconds for case 1. 

Dimensions as in figure 1. Areas in which Ri < 1 are 
indicated by a dashed line. Areas in which Ri < 1/4 
are indicated by the symbol "A". 
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Figure 9. 


Streamfunction field at 3750 seconds 
shear, constant stability). Ax = 2000 
horizontal extent displayed = 128 km, 
displayed = 20 km. 
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Figure 1 0. 


Streamfunction field at 7500 seconds for case 2 , 
as in figure 9, 
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Figure 1 1 . 


Streamfunction field at 1200 seconds for 

1 /2 

shear, constant stability, Ri^ ' =3,3). 
500 m, horizontal extent displayed = 32 
displayed = 7.5 km. 
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Figure 12. Streamfunction field at 4500 seconds for case 4 (exponential 

1 /2 

shear, constant stability, Ri^ =6.0). Ax = 750 m, Az = 
500 m, horizontal extent displayed = 48 km, vertical extent 
displayed = 7.5 km. 
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Figure 13. Streamfunction field at 3000 seconds for case 5 (nonlinear, 

2 

constant pu ) . Ax = 1000 m, Az = 1000 m, horizontal extent 
displayed = 64 km, vertical extent displayed = 10 km. 
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Figure 14, Vertical velocity field at 3000 seconds for case 5, 
Dimensions as in figure 13, 
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Figure 16. 


Density 


field at 1500 seconds 


for case 6. 


figure 15, 


Dimensions as in 
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Figure 17, 


Streamfunction field at 600 seconds for case 7 (constant 
velocity, constant stability, large obstacle, k^h=1.95). 

Ax = 625 m, Az = 625 m, horizontal extent displayed =40 km, 
vertical extent displayed = 14,38 km. 
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Figure 18. 


Streamfunction field 


at 800 seconds 


for case 7 . 


as i n figure 17 . 
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Figure 19. 
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field at 1000 seconds 
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as in figure 17. 
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Figure 20. 


Density field at 1000 seconds for case 7 , 
figure 17, 


Dimensions as in 
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Richardson number field at 1 
Dimensions as in figure 15, 
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Figure 22. Richardson number field at 1000 seconds for case 7. 

Dimensions as in figure 17. Areas in which Ri < 1 are 
indicated by a dashed line. Areas in which Ri < 1/4 
are indicated by the symbol "A". 
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Figure 23. Cross section of the potential temperature (in K) along 
an east-west line through Boulder on 11 January 1972. 
Analysis above the heavy dashed line is from the Sabreliner 
data, taken between 1700 and 2000 MST, and analysis below 
this line is primarily from the Queen Air data, taken 
between 1330 and 1500 MST. Flight tracks are indicated by 
the light dashed lines (from Lilly and Zipser, 1972). 
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Figure 24, Cross section of horizontal wind velocity (in m/sec) along 
an east-west line through Boulder on 11 January 1972. This 
analysis was derived from Sabreliner data only. The analysis 
below 500 mb was partially obtained from vertical integration 
of the continuity equation, assuming two-dimensional, steady- 
state flow. Crosses indicate turbulent portions along the 
flight track (from Lilly and Zipser, 1972). 
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Figure 25. 
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Figure 26. 


Horizontal velocity field at 4250 seconds for case 8 


Dimensions as 


in figure 25. 
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